model.py 7.47 KB
 Tomas Pettersson committed Nov 24, 2017 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ``````# Particles are distributed at the nodes of a grid formed by equilateral # triangles. The distance between the particles, ie the length of the sides of # the triangles,is determined by the requirement that the area one particle # represents times the number of particles should equal the area of the polygon # determined by the outlet vertices. The "area" of one particle corresponds to # twice the area of one equilateral triangle, which equals # sqrt(3)/4*ParticleDistance**2. # A grid of equilateral triangles forms regular hexagons. The algorithm to # distribute particles is to determine one starting point, then move around # this following the sides of increasingly larger regular hexagons. The first # hexagon around the starting point has 6 vertices, the next has 12, the next # has 18, and so on. # The number of nodes are # 1+6*1+6*2+6*3+...+6*n = 1+6*(1+2+3...+n) = 1+6*(n/2*(n+1)) # where n is the number of hexagons (n is also the number of ParticleDistances # from the center to a node on the outermost hexagon) import math import json import random `````` Tomas Pettersson committed Dec 07, 2017 22 23 24 ``````from output import Output from output import Particletrack from output import Cloudtrack `````` Tomas Pettersson committed Nov 24, 2017 25 26 27 28 29 30 31 32 33 34 35 36 37 38 ``````from shapely import geometry class Model(object): def __init__(self): return def __call__(self): return self def centerPoints(self, newLatLng, centroid, points): centered = [] center = [centroid.x,centroid.y] for point in points: latlng = [point.x, point.y] `````` Tomas Pettersson committed Nov 29, 2017 39 `````` centered.append(self.centerPoint(newLatLng, center, latlng)) `````` Tomas Pettersson committed Nov 24, 2017 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 `````` return centered def centerPoint(self, newCenter, oldCenter, latlng): if (latlng[0] > oldCenter[0]): lat = abs(latlng[0] - oldCenter[0]) + newCenter[0] else: lat = newCenter[0] - abs(latlng[0] - oldCenter[0]) if (latlng[1] > oldCenter[1]): lng = abs(latlng[1] - oldCenter[1]) + newCenter[1] else: lng = newCenter[1] - abs(latlng[1] - oldCenter[1]) return [lat, lng] def createHexagon(self, size, center, lvl): sizeLvl = size * lvl coords = [] nodes = [] for i in range(6): angle = 2 * math.pi / 6 * (i + 0.5) lat_i = center.x + sizeLvl * math.cos(angle) lng_i = center.y + sizeLvl * math.sin(angle) * 1.5 coords.append(geometry.Point(lat_i, lng_i)) for i, coord in enumerate(coords): nodes.append(coord) if lvl > 1: index = 0 if (i + 1) < len(coords): index = i + 1 latDiff = (coords[index].x - coord.x) / lvl lngDiff = (coords[index].y - coord.y) / lvl j = 1 while (j < lvl): nodes.append(geometry.Point(coord.x + (latDiff * j), coord.y + (lngDiff * j))) j += 1 if (lvl == 1): nodes.append(center) return { 'center': center, 'coordinates': coords, 'nodes': nodes } `````` Tomas Pettersson committed Dec 07, 2017 88 `````` def createOutlet(self,geom): `````` Tomas Pettersson committed Dec 12, 2017 89 `````` nrOfParticles = 300 `````` Tomas Pettersson committed Dec 07, 2017 90 `````` center = geom.centroid `````` Tomas Pettersson committed Nov 24, 2017 91 92 93 94 `````` points = [] depth = [] properties = {} properties['depth'] = depth `````` Tomas Pettersson committed Dec 07, 2017 95 `````` if (geom.geom_type == 'Point'): `````` Tomas Pettersson committed Dec 13, 2017 96 `````` geom = geom.buffer(0.001) `````` Tomas Pettersson committed Dec 07, 2017 97 98 99 `````` if (geom.geom_type == 'LineString'): geom = geom.buffer(0.01, 20) pArea = (geom.area)/3 `````` Tomas Pettersson committed Nov 24, 2017 100 101 102 103 104 105 106 107 108 `````` pDist = math.sqrt((4*pArea)/(nrOfParticles*math.sqrt(3))) counter = 0 lvl = 1 while (counter < nrOfParticles): hexagon = self.createHexagon(pDist, center, lvl) nodes = hexagon['nodes'] j = 0 while (j < len(nodes)): point = nodes[j] `````` Tomas Pettersson committed Dec 07, 2017 109 `````` if geom.contains(point): `````` Tomas Pettersson committed Nov 24, 2017 110 111 112 113 114 115 116 117 118 119 `````` points.append(point) depth.append(0.0) counter += 1 if (counter == nrOfParticles): break j += 1 lvl += 1 if (lvl > 100): break `````` Tomas Pettersson committed Nov 29, 2017 120 `````` print("Outlet particles count: "+str(counter)) `````` Tomas Pettersson committed Nov 24, 2017 121 `````` `````` Tomas Pettersson committed Dec 12, 2017 122 `````` return geometry.MultiPoint(points),pDist `````` Tomas Pettersson committed Nov 24, 2017 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 `````` def createFeatureCollection(self): featurecollection = {} featurecollection["type"] = "FeatureCollection" featurecollection["features"] = [] return featurecollection def createFeature(self, geom, properties = {}): featurecollection = {} featurecollection["type"] = "Feature" featurecollection["properties"] = properties featurecollection["geometry"] = geometry.mapping(geom) return featurecollection def onlandPoints(self, points, strtree): result = list(points) `````` Tomas Pettersson committed Nov 29, 2017 140 `````` print("On land calculated") `````` Tomas Pettersson committed Nov 24, 2017 141 142 `````` return result `````` Tomas Pettersson committed Dec 12, 2017 143 144 `````` def calculateRadius(self, pDist, points): radius = [pDist] * len(points) `````` Tomas Pettersson committed Nov 24, 2017 145 `````` `````` Tomas Pettersson committed Dec 12, 2017 146 147 148 149 `````` for i in range(len(points)): radius[i] *= random.uniform(0.05, 3) return radius `````` Tomas Pettersson committed Nov 24, 2017 150 151 `````` `````` Tomas Pettersson committed Dec 12, 2017 152 153 `````` def displacePoints(self, pDist, points): `````` Tomas Pettersson committed Nov 24, 2017 154 155 156 157 158 `````` p = 0 xdisp = [0] * len(points) ydisp = [0] * len(points) minoverlap = [float('Inf')] * len(points) `````` Tomas Pettersson committed Dec 12, 2017 159 `````` radius = self.calculateRadius(pDist, points) `````` Tomas Pettersson committed Nov 24, 2017 160 `````` while (p < len(points)): `````` Tomas Pettersson committed Dec 12, 2017 161 `````` pn = 0 `````` Tomas Pettersson committed Nov 24, 2017 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 `````` while (pn < len(points)): dx = points[p][0] - points[pn][0] dy = points[p][1] - points[pn][1] centerdist = math.sqrt(dx*dx + dy*dy) overlap = radius[p] + radius[pn] - centerdist if overlap > 0: minoverlap[p] = min(minoverlap[p], overlap) minoverlap[pn] = min(minoverlap[pn], overlap) if centerdist > 0: cosalfa = (points[pn][0] - points[p][0]) / centerdist sinalfa = (points[pn][1] - points[p][1]) / centerdist else: rand = random.uniform(0, 1) cosalfa = math.cos(rand * math.pi*math.pi) sinalfa = math.sin(rand * math.pi*math.pi) xdisp[p] = xdisp[p] - 0.5 * cosalfa * overlap ydisp[p] = ydisp[p] - 0.5 * sinalfa * overlap xdisp[pn] = xdisp[pn] + 0.5 * cosalfa * overlap ydisp[pn] = ydisp[pn] + 0.5 * sinalfa * overlap pn += 1 disp = math.sqrt(xdisp[p]*xdisp[p]+ydisp[p]*ydisp[p]) maxdisp = min(radius[p], 0.5 * minoverlap[p]) if (disp > maxdisp): xdisp[p] = xdisp[p] * maxdisp / disp ydisp[p] = ydisp[p] * maxdisp / disp p += 1 result = list(points) for i, dx in enumerate(xdisp): result[i][0] += dx for i, dy in enumerate(ydisp): result[i][1] += dy return result def createProperties(self, step, time, centroid, level, category): properties = {}; properties['nStep'] = step properties['time'] = time `````` Tomas Pettersson committed Nov 29, 2017 206 207 `````` properties['meanLat'] = centroid.y properties['meanLon'] = centroid.x `````` Tomas Pettersson committed Nov 24, 2017 208 209 210 211 212 `````` properties['category'] = category properties['level'] = level return properties `````` Tomas Pettersson committed Dec 12, 2017 213 `````` def createOutput(self, starttime, pDist, multipoint, exercisefeature, strtree): `````` Tomas Pettersson committed Nov 24, 2017 214 215 216 217 `````` featurecollection = self.createFeatureCollection() features = featurecollection['features'] time = exercisefeature['properties']['time'] `````` Tomas Pettersson committed Dec 06, 2017 218 `````` properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms)) `````` Tomas Pettersson committed Nov 24, 2017 219 220 221 222 223 `````` features.append(self.createFeature(multipoint, properties)) mp = multipoint linestring = geometry.shape(exercisefeature['geometry']) for i, coord in enumerate(linestring.coords): points = self.centerPoints(coord, mp.centroid, mp.geoms) `````` Tomas Pettersson committed Dec 12, 2017 224 `````` displacedpoints = self.displacePoints(pDist, points) `````` Tomas Pettersson committed Nov 24, 2017 225 `````` # onlandpoints = self.onlandPoints(displacedpoints,strtree) `````` Tomas Pettersson committed Dec 12, 2017 226 `````` mp = geometry.MultiPoint(displacedpoints) `````` Tomas Pettersson committed Nov 24, 2017 227 228 `````` level = [0] * len(mp.geoms) category = [2] * len(mp.geoms) `````` Tomas Pettersson committed Dec 06, 2017 229 `````` properties = self.createProperties((i+2), time[i], mp.centroid, level, category) `````` Tomas Pettersson committed Nov 24, 2017 230 `````` features.append(self.createFeature(mp, properties)) `````` Tomas Pettersson committed Dec 07, 2017 231 232 233 234 `````` if (i % 10 == 0): Cloudtrack.write(featurecollection) Particletrack.write(featurecollection) `````` Tomas Pettersson committed Nov 24, 2017 235 `` return featurecollection``