# 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 from output import Output from output import Particletrack from output import Cloudtrack 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] 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, points, strtree): result = list(points) print("On land calculated") return result 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'] time = 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 linestring = geometry.shape(exercisefeature['geometry']) for i, coord in enumerate(linestring.coords): points = self.centerPoints(coord, mp.centroid, mp.geoms) displacedpoints = self.displacePoints(pDist, points) # onlandpoints = self.onlandPoints(displacedpoints,strtree) mp = geometry.MultiPoint(displacedpoints) level = [0] * len(mp.geoms) category = [2] * len(mp.geoms) properties = self.createProperties((i+2), time[i], mp.centroid, level, category) features.append(self.createFeature(mp, properties)) if (i % 10 == 0): Cloudtrack.write(featurecollection) Particletrack.write(featurecollection) return featurecollection