# 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 import time from output import Output from output import Particletrack from output import Cloudtrack from shapely import geometry from shapely.geometry import LineString from shapely.strtree import STRtree from shapely.geometry import MultiLineString 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] centered.append(self.centerPoint(newLatLng, center, latlng)) 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 } def createOutlet(self,geom): nrOfParticles = 300 center = geom.centroid points = [] depth = [] properties = {} properties['depth'] = depth if (geom.geom_type == 'Point'): geom = geom.buffer(0.001) if (geom.geom_type == 'LineString'): geom = geom.buffer(0.01, 20) pArea = (geom.area)/3 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] if geom.contains(point): points.append(point) depth.append(0.0) counter += 1 if (counter == nrOfParticles): break j += 1 lvl += 1 if (lvl > 100): break print("Outlet particles count: "+str(counter)) return geometry.MultiPoint(points),pDist 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, 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 def calculateRadius(self, pDist, points): radius = [pDist] * len(points) for i in range(len(points)): radius[i] *= random.uniform(0.05, 3) return radius def displacePoints(self, pDist, points): p = 0 xdisp = [0] * len(points) ydisp = [0] * len(points) minoverlap = [float('Inf')] * len(points) radius = self.calculateRadius(pDist, points) while (p < len(points)): pn = 0 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 properties['meanLat'] = centroid.y properties['meanLon'] = centroid.x properties['category'] = category properties['level'] = level return properties def createOutput(self, starttime, pDist, multipoint, exercisefeature, strtree): featurecollection = self.createFeatureCollection() features = featurecollection['features'] timearray = exercisefeature['properties']['time'] properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms)) features.append(self.createFeature(multipoint, properties)) mp = multipoint activepoints = [] deactivepoints = [] linestring = geometry.shape(exercisefeature['geometry']) latlngpoints = [[latlng[1], latlng[0]] for latlng in linestring.coords] print('before matches') matches = strtree.query(LineString(latlngpoints)) print('matches: '+str(matches)) for i, coord in enumerate(linestring.coords): start = time.time() 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))) end = time.time() print('onland: '+str(end-start)) mp = geometry.MultiPoint(activepoints+deactivepoints) level = [0] * len(mp.geoms) category = [2] * len(mp.geoms) properties = self.createProperties((i+2), timearray[i], mp.centroid, level, category) features.append(self.createFeature(mp, properties)) if (i % 10 == 0): Cloudtrack.write(featurecollection) Particletrack.write(featurecollection) mp = geometry.MultiPoint(activepoints) print('deactivepoints: '+str(len(deactivepoints))) return featurecollection