model.py 7.38 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 Nov 24, 2017 89 `````` nrOfParticles = 500 `````` 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 96 97 98 99 `````` if (geom.geom_type == 'Point'): geom = geom.buffer(0.0001) 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 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 `````` return geometry.MultiPoint(points) 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 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 `````` return result def displacePoints(self, points): p = 0 xdisp = [0] * len(points) ydisp = [0] * len(points) minoverlap = [float('Inf')] * len(points) radius = [0] * len(points) for i in range(len(points)): radius[i] = 0.035 * random.uniform(0.1, 1) while (p < len(points)): pn = p + 1 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 201 202 `````` properties['meanLat'] = centroid.y properties['meanLon'] = centroid.x `````` Tomas Pettersson committed Nov 24, 2017 203 204 205 206 207 `````` properties['category'] = category properties['level'] = level return properties `````` Tomas Pettersson committed Dec 06, 2017 208 `````` def createOutput(self, starttime, multipoint, exercisefeature, strtree): `````` Tomas Pettersson committed Nov 24, 2017 209 210 211 212 `````` featurecollection = self.createFeatureCollection() features = featurecollection['features'] time = exercisefeature['properties']['time'] `````` Tomas Pettersson committed Dec 06, 2017 213 `````` properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms)) `````` Tomas Pettersson committed Nov 24, 2017 214 215 216 217 218 `````` 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 11, 2017 219 `````` # displacedpoints = self.displacePoints(points) `````` Tomas Pettersson committed Nov 24, 2017 220 `````` # onlandpoints = self.onlandPoints(displacedpoints,strtree) `````` Tomas Pettersson committed Dec 11, 2017 221 `````` mp = geometry.MultiPoint(points) `````` Tomas Pettersson committed Nov 24, 2017 222 223 `````` level = [0] * len(mp.geoms) category = [2] * len(mp.geoms) `````` Tomas Pettersson committed Dec 06, 2017 224 `````` properties = self.createProperties((i+2), time[i], mp.centroid, level, category) `````` Tomas Pettersson committed Nov 24, 2017 225 `````` features.append(self.createFeature(mp, properties)) `````` Tomas Pettersson committed Dec 07, 2017 226 227 228 229 230 `````` if (i % 10 == 0): Cloudtrack.write(featurecollection) Particletrack.write(featurecollection) Output.write(featurecollection) `````` Tomas Pettersson committed Nov 24, 2017 231 `` return featurecollection``