model.py 9.59 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 13, 2017 22 ``````import time `````` Tomas Pettersson committed Dec 07, 2017 23 24 25 ``````from output import Output from output import Particletrack from output import Cloudtrack `````` Tomas Pettersson committed Nov 24, 2017 26 ``````from shapely import geometry `````` Tomas Pettersson committed Dec 13, 2017 27 28 29 30 ``````from shapely.geometry import LineString from shapely.strtree import STRtree from shapely.geometry import MultiLineString `````` Tomas Pettersson committed Nov 24, 2017 31 32 33 34 35 36 37 38 39 40 41 42 43 `````` 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 44 `````` centered.append(self.centerPoint(newLatLng, center, latlng)) `````` Tomas Pettersson committed Nov 24, 2017 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 88 89 90 91 92 `````` 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 93 `````` def createOutlet(self,geom): `````` Tomas Pettersson committed Dec 14, 2017 94 `````` nrOfParticles = 300 `````` Tomas Pettersson committed Dec 07, 2017 95 `````` center = geom.centroid `````` Tomas Pettersson committed Nov 24, 2017 96 97 98 99 `````` points = [] depth = [] properties = {} properties['depth'] = depth `````` Tomas Pettersson committed Dec 07, 2017 100 `````` if (geom.geom_type == 'Point'): `````` Tomas Pettersson committed Dec 13, 2017 101 `````` geom = geom.buffer(0.001) `````` Tomas Pettersson committed Dec 07, 2017 102 103 104 `````` if (geom.geom_type == 'LineString'): geom = geom.buffer(0.01, 20) pArea = (geom.area)/3 `````` Tomas Pettersson committed Nov 24, 2017 105 106 107 108 109 110 111 112 113 `````` 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 114 `````` if geom.contains(point): `````` Tomas Pettersson committed Nov 24, 2017 115 116 117 118 119 120 121 122 123 124 `````` 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 125 `````` print("Outlet particles count: "+str(counter)) `````` Tomas Pettersson committed Nov 24, 2017 126 `````` `````` Tomas Pettersson committed Dec 12, 2017 127 `````` return geometry.MultiPoint(points),pDist `````` Tomas Pettersson committed Nov 24, 2017 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 `````` 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 `````` Tomas Pettersson committed Dec 15, 2017 143 `````` def shorePoints(self, pointsdata, beforepoints, nowpoints, activeindex, matches): `````` Tomas Pettersson committed Dec 13, 2017 144 145 146 147 148 149 150 151 `````` for i in range(len(beforepoints)): queryline = LineString([[beforepoints[i].y,beforepoints[i].x],[nowpoints[i][1],nowpoints[i][0]]]) intersection = None if (len(matches) > 0): for match in matches: if match.intersects(queryline): intersection = match.intersection(queryline) if (intersection is not None): `````` Tomas Pettersson committed Dec 15, 2017 152 `````` pointsdata[activeindex[i]][1] = 7 # on shore category `````` Tomas Pettersson committed Dec 14, 2017 153 `````` if (intersection.geom_type is 'Point'): `````` Tomas Pettersson committed Dec 15, 2017 154 155 `````` pointsdata[activeindex[i]][0] = [intersection.y,intersection.x] # deactive.append([intersection.y,intersection.x]) `````` Tomas Pettersson committed Dec 14, 2017 156 `````` if (intersection.geom_type is 'MultiPoint'): `````` Tomas Pettersson committed Dec 15, 2017 157 158 `````` pointsdata[activeindex[i]][0] = [intersection.geoms[0].y,intersection.geoms[0].x] # deactive.append([intersection.geoms[0].y,intersection.geoms[0].x]) `````` Tomas Pettersson committed Dec 13, 2017 159 `````` `````` Tomas Pettersson committed Dec 15, 2017 160 `````` return pointsdata `````` Tomas Pettersson committed Dec 13, 2017 161 `````` `````` Tomas Pettersson committed Nov 24, 2017 162 `````` `````` Tomas Pettersson committed Dec 12, 2017 163 164 `````` def calculateRadius(self, pDist, points): radius = [pDist] * len(points) `````` Tomas Pettersson committed Nov 24, 2017 165 `````` `````` Tomas Pettersson committed Dec 12, 2017 166 167 168 169 `````` for i in range(len(points)): radius[i] *= random.uniform(0.05, 3) return radius `````` Tomas Pettersson committed Nov 24, 2017 170 `````` `````` Tomas Pettersson committed Dec 12, 2017 171 `````` def displacePoints(self, pDist, points): `````` Tomas Pettersson committed Nov 24, 2017 172 173 174 175 176 `````` p = 0 xdisp = [0] * len(points) ydisp = [0] * len(points) minoverlap = [float('Inf')] * len(points) `````` Tomas Pettersson committed Dec 12, 2017 177 `````` radius = self.calculateRadius(pDist, points) `````` Tomas Pettersson committed Nov 24, 2017 178 `````` while (p < len(points)): `````` Tomas Pettersson committed Dec 12, 2017 179 `````` pn = 0 `````` Tomas Pettersson committed Nov 24, 2017 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 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 `````` 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 224 225 `````` properties['meanLat'] = centroid.y properties['meanLon'] = centroid.x `````` Tomas Pettersson committed Nov 24, 2017 226 227 228 229 `````` properties['category'] = category properties['level'] = level return properties `````` Tomas Pettersson committed Dec 15, 2017 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 `````` def updatepointarrays(self, pointsdata): activepoints = [] deactivepoints = [] activeindex = [] allpoints = [] category = [] for i,pointdata in enumerate(pointsdata): allpoints.append(pointdata[0]) category.append(pointdata[1]) if (category[i] is 2): # active activepoints.append(pointdata[0]) activeindex.append(i) else: # not active deactivepoints.append(pointdata[0]) return allpoints, activepoints, activeindex, deactivepoints, category `````` Tomas Pettersson committed Nov 24, 2017 245 `````` `````` Tomas Pettersson committed Dec 12, 2017 246 `````` def createOutput(self, starttime, pDist, multipoint, exercisefeature, strtree): `````` Tomas Pettersson committed Nov 24, 2017 247 248 249 `````` featurecollection = self.createFeatureCollection() features = featurecollection['features'] `````` Tomas Pettersson committed Dec 13, 2017 250 `````` timearray = exercisefeature['properties']['time'] `````` Tomas Pettersson committed Dec 06, 2017 251 `````` properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms)) `````` Tomas Pettersson committed Nov 24, 2017 252 `````` features.append(self.createFeature(multipoint, properties)) `````` Tomas Pettersson committed Dec 15, 2017 253 254 255 256 257 `````` pointsdata = () for i,point in enumerate(multipoint.geoms): pointsdata += ([[point.x, point.y], 2],) # default category 2 = active in water allpoints, activepoints, activeindex, deactivepoints, category = self.updatepointarrays(pointsdata) `````` Tomas Pettersson committed Nov 24, 2017 258 `````` linestring = geometry.shape(exercisefeature['geometry']) `````` Tomas Pettersson committed Dec 13, 2017 259 `````` latlngpoints = [[latlng[1], latlng[0]] for latlng in linestring.coords] `````` Tomas Pettersson committed Dec 14, 2017 260 `````` matches = strtree.query(LineString(latlngpoints)) `````` Tomas Pettersson committed Nov 24, 2017 261 `````` for i, coord in enumerate(linestring.coords): `````` Tomas Pettersson committed Dec 14, 2017 262 263 `````` if (len(activepoints) > 0): mp = geometry.MultiPoint(activepoints) `````` Tomas Pettersson committed Dec 15, 2017 264 265 266 267 268 `````` centeredpoints = self.centerPoints(coord, mp.centroid, mp.geoms) displacedpoints = self.displacePoints(pDist, centeredpoints) pointsdata = self.shorePoints(pointsdata, mp.geoms, displacedpoints, activeindex, matches) allpoints, activepoints, activeindex, deactivepoints, category = self.updatepointarrays(pointsdata) mp = geometry.MultiPoint(allpoints) `````` Tomas Pettersson committed Nov 24, 2017 269 `````` level = [0] * len(mp.geoms) `````` Tomas Pettersson committed Dec 13, 2017 270 `````` properties = self.createProperties((i+2), timearray[i], mp.centroid, level, category) `````` Tomas Pettersson committed Nov 24, 2017 271 `````` features.append(self.createFeature(mp, properties)) `````` Tomas Pettersson committed Dec 14, 2017 272 273 274 `````` if (i % 10 == 0): Cloudtrack.write(featurecollection) Particletrack.write(featurecollection) `````` Tomas Pettersson committed Dec 14, 2017 275 `````` `````` Tomas Pettersson committed Dec 14, 2017 276 `````` print('activepoints: '+str(len(activepoints))) `````` Tomas Pettersson committed Dec 14, 2017 277 `````` print('deactivepoints: '+str(len(deactivepoints))) `````` Tomas Pettersson committed Dec 07, 2017 278 `````` `````` Tomas Pettersson committed Nov 24, 2017 279 `` return featurecollection``