model.py 10.2 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 13, 2017 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 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 `````` def onlandPoints(self, beforepoints, nowpoints, matches): active = [] deactive = [] 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): deactive.append([intersection.coords[0][0],intersection.coords[0][1]]) else: active.append(nowpoints[i]) intersectionend = time.time() return active, deactive def onlandPoints2(self, beforepoints, nowpoints, matches): active = [] deactive = [] multilines = [] intersectionstart = time.time() for i in range(len(beforepoints)): multilines.append(LineString([[beforepoints[i].y,beforepoints[i].x],[nowpoints[i][1],nowpoints[i][0]]])) intercepts = [match.intersection(MultiLineString(multilines)) for match in matches] # if (len(matches) == 3): # print(len(matches[2].coords)) # print(matches[2].intersection(query)) # intercepts = [match.intersection(query) for match in matches] # if (len(intercepts)): # print(intercepts) # for i in range(len(multilines)): # intercepts = [multilines[i] for match in matches if match.intersects(multilines[i])] # if (len(intercepts)): # print(intercepts) active = list(nowpoints) # if (len(matches) > 0): # print(len(matches)) # for i in range(len(multilines)): # intersection = None # for match in matches: # intersection = match.intersection(multilines[i]) # if (not intersection.is_empty): # break # # if match.intersects(multilines[i]): # print(intersection) # if (intersection.is_empty): # active.append(nowpoints[i]) # else: # deactive.append([intersection.y,intersection.x]) intersectionend = time.time() print('intersection: '+str(intersectionend-intersectionstart)) return active, deactive `````` Tomas Pettersson committed Nov 24, 2017 205 `````` `````` Tomas Pettersson committed Dec 12, 2017 206 207 `````` def calculateRadius(self, pDist, points): radius = [pDist] * len(points) `````` Tomas Pettersson committed Nov 24, 2017 208 `````` `````` Tomas Pettersson committed Dec 12, 2017 209 210 211 212 `````` for i in range(len(points)): radius[i] *= random.uniform(0.05, 3) return radius `````` Tomas Pettersson committed Nov 24, 2017 213 214 `````` `````` Tomas Pettersson committed Dec 12, 2017 215 216 `````` def displacePoints(self, pDist, points): `````` Tomas Pettersson committed Nov 24, 2017 217 218 219 220 221 `````` p = 0 xdisp = [0] * len(points) ydisp = [0] * len(points) minoverlap = [float('Inf')] * len(points) `````` Tomas Pettersson committed Dec 12, 2017 222 `````` radius = self.calculateRadius(pDist, points) `````` Tomas Pettersson committed Nov 24, 2017 223 `````` while (p < len(points)): `````` Tomas Pettersson committed Dec 12, 2017 224 `````` pn = 0 `````` Tomas Pettersson committed Nov 24, 2017 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 `````` 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 269 270 `````` properties['meanLat'] = centroid.y properties['meanLon'] = centroid.x `````` Tomas Pettersson committed Nov 24, 2017 271 272 273 274 275 `````` properties['category'] = category properties['level'] = level return properties `````` Tomas Pettersson committed Dec 12, 2017 276 `````` def createOutput(self, starttime, pDist, multipoint, exercisefeature, strtree): `````` Tomas Pettersson committed Nov 24, 2017 277 278 279 `````` featurecollection = self.createFeatureCollection() features = featurecollection['features'] `````` Tomas Pettersson committed Dec 13, 2017 280 `````` timearray = exercisefeature['properties']['time'] `````` Tomas Pettersson committed Dec 06, 2017 281 `````` properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms)) `````` Tomas Pettersson committed Nov 24, 2017 282 283 `````` features.append(self.createFeature(multipoint, properties)) mp = multipoint `````` Tomas Pettersson committed Dec 14, 2017 284 `````` activepoints = [] `````` Tomas Pettersson committed Dec 13, 2017 285 `````` deactivepoints = [] `````` Tomas Pettersson committed Nov 24, 2017 286 `````` linestring = geometry.shape(exercisefeature['geometry']) `````` Tomas Pettersson committed Dec 13, 2017 287 `````` latlngpoints = [[latlng[1], latlng[0]] for latlng in linestring.coords] `````` Tomas Pettersson committed Dec 14, 2017 288 `````` print('before matches') `````` Tomas Pettersson committed Dec 13, 2017 289 `````` matches = strtree.query(LineString(latlngpoints)) `````` Tomas Pettersson committed Dec 14, 2017 290 `````` print('matches: '+str(matches)) `````` Tomas Pettersson committed Nov 24, 2017 291 `````` for i, coord in enumerate(linestring.coords): `````` Tomas Pettersson committed Dec 14, 2017 292 `````` start = time.time() `````` Tomas Pettersson committed Dec 14, 2017 293 294 295 296 297 298 299 `````` if (len(mp.geoms) > 0): points = self.centerPoints(coord, mp.centroid, mp.geoms) displacedpoints = self.displacePoints(pDist, points) activepoints, onlandpoints = self.onlandPoints(mp.geoms, displacedpoints, matches) deactivepoints += onlandpoints # print('active: '+str(len(activepoints))) # print('onland: '+str(len(onlandpoints))) `````` Tomas Pettersson committed Dec 14, 2017 300 301 `````` end = time.time() print('onland: '+str(end-start)) `````` Tomas Pettersson committed Dec 13, 2017 302 `````` mp = geometry.MultiPoint(activepoints+deactivepoints) `````` Tomas Pettersson committed Nov 24, 2017 303 304 `````` level = [0] * len(mp.geoms) category = [2] * len(mp.geoms) `````` Tomas Pettersson committed Dec 13, 2017 305 `````` properties = self.createProperties((i+2), timearray[i], mp.centroid, level, category) `````` Tomas Pettersson committed Nov 24, 2017 306 `````` features.append(self.createFeature(mp, properties)) `````` Tomas Pettersson committed Dec 14, 2017 307 308 309 `````` if (i % 10 == 0): Cloudtrack.write(featurecollection) Particletrack.write(featurecollection) `````` Tomas Pettersson committed Dec 14, 2017 310 311 312 `````` mp = geometry.MultiPoint(activepoints) print('deactivepoints: '+str(len(deactivepoints))) `````` Tomas Pettersson committed Dec 07, 2017 313 `````` `````` Tomas Pettersson committed Nov 24, 2017 314 `` return featurecollection``