Commit 57250af5 authored by Tomas Pettersson's avatar Tomas Pettersson 🏸

onland fixes

parent cba5dd25
2017-12-12 12:00 19.7563 57.5310 0.0000 2017-12-13 08:00 27.0119 60.0031 0.0000
2017-12-12 12:15 19.7699 57.5309 0.0000 2017-12-13 08:15 27.0112 60.0034 0.0000
2017-12-12 12:30 19.7780 57.5261 0.0000
2017-12-12 12:45 19.7903 57.5208 0.0000
2017-12-12 13:00 19.7997 57.5166 0.0000
2017-12-12 13:15 19.8109 57.5116 0.0000
2017-12-12 13:30 19.8182 57.5084 0.0000
2017-12-12 13:45 19.8274 57.5047 0.0000
2017-12-12 14:00 19.8399 57.5006 0.0000
2017-12-12 14:15 19.8564 57.4959 0.0000
2017-12-12 14:30 19.8688 57.4930 0.0000
2017-12-12 14:45 19.8823 57.4922 0.0000
2017-12-12 15:00 19.8920 57.4921 0.0000
2017-12-12 15:15 19.9030 57.4931 0.0000
2017-12-12 15:30 19.9113 57.4953 0.0000
2017-12-12 15:45 19.9194 57.4988 0.0000
2017-12-12 16:00 19.9291 57.5025 0.0000
2017-12-12 16:15 19.9389 57.5078 0.0000
2017-12-12 16:30 19.9471 57.5128 0.0000
2017-12-12 16:45 19.9539 57.5174 0.0000
2017-12-12 17:00 19.9607 57.5209 0.0000
2017-12-12 17:15 19.9675 57.5254 0.0000
2017-12-12 17:30 19.9730 57.5290 0.0000
2017-12-12 17:45 19.9800 57.5319 0.0000
2017-12-12 18:00 19.9909 57.5342 0.0000
2017-12-12 18:15 20.0048 57.5358 0.0000
2017-12-12 18:30 20.0143 57.5372 0.0000
2017-12-12 18:45 20.0239 57.5380 0.0000
2017-12-12 19:00 20.0377 57.5379 0.0000
2017-12-12 19:15 20.0514 57.5379 0.0000
2017-12-12 19:30 20.0611 57.5380 0.0000
2017-12-12 19:45 20.0775 57.5366 0.0000
2017-12-12 20:00 20.0871 57.5342 0.0000
2017-12-12 20:15 20.0981 57.5299 0.0000
2017-12-12 20:30 20.1063 57.5254 0.0000
2017-12-12 20:45 20.1118 57.5209 0.0000
2017-12-12 21:00 20.1160 57.5159 0.0000
2017-12-12 21:15 20.1187 57.5115 0.0000
2017-12-12 21:30 20.1201 57.5070 0.0000
2017-12-12 21:45 20.1201 57.4981 0.0000
2017-12-12 22:00 20.1200 57.4922 0.0000
2017-12-12 22:15 20.1228 57.4870 0.0000
2017-12-12 22:30 20.1229 57.4819 0.0000
2017-12-12 22:45 20.1173 57.4768 0.0000
2017-12-12 23:00 20.1090 57.4730 0.0000
2017-12-12 23:15 20.1022 57.4692 0.0000
2017-12-12 23:30 20.0967 57.4649 0.0000
2017-12-12 23:45 20.0912 57.4596 0.0000
2017-12-13 00:00 20.0871 57.4530 0.0000
...@@ -16,20 +16,35 @@ from shapely.geometry import LineString ...@@ -16,20 +16,35 @@ from shapely.geometry import LineString
from shapely.strtree import STRtree from shapely.strtree import STRtree
def storetree(): def storetree():
with open('../coastline_HELCOM_82214.dat',"r") as coastlinefile: try:
lines = coastlinefile.readlines() with open('../coastline_HELCOM_82214.dat',"r") as coastlinefile:
section = [] lines = coastlinefile.readlines()
sections = [] section = []
for line in lines: sections = []
if '-999999' in line:
if len(section) > 0: for line in lines:
if len(section) == 10:
sections.append(LineString(section)) sections.append(LineString(section))
section = [] section = []
continue if '-999999' in line:
split = line.split(' ') if len(section) > 1: #this will skip some points where section will end with only one coord. Possibly create a line from previous section on ending point
section.append([float(split[2]),float(split[1])]) sections.append(LineString(section))
with open('resources/strtree.pickle', 'wb') as handle: section = []
pickle.dump(sections, handle, protocol=pickle.HIGHEST_PROTOCOL) continue
split = line.split(' ')
section.append([float(split[2]),float(split[1])])
print('Number of lines: '+str(len(sections)))
with open('resources/strtree.pickle', 'wb') as handle:
pickle.dump(sections, handle, protocol=pickle.HIGHEST_PROTOCOL)
# verify
tree = STRtree(sections)
queryline = LineString([[57.5,17.5],[57.5,18.5]])
matches = tree.query(queryline);
print('Number of matches: '+str(len(matches)))
except Exception as error:
print_exc()
sys.exit(1)
def read_properties(filepath, sep='=', comment_char='#'): def read_properties(filepath, sep='=', comment_char='#'):
...@@ -85,8 +100,9 @@ def run(): ...@@ -85,8 +100,9 @@ def run():
outputfeaturecollection['properties']['status'] = 'COMPLETE' outputfeaturecollection['properties']['status'] = 'COMPLETE'
outputfeaturecollection['properties']['simulation'] = outletproperties['simulation'] outputfeaturecollection['properties']['simulation'] = outletproperties['simulation']
Cloudtrack.write(outputfeaturecollection) #Cloudtrack.write(outputfeaturecollection)
Particletrack.write(outputfeaturecollection) #Particletrack.write(outputfeaturecollection)
Output.write(outputfeaturecollection)
print("COMPLETE") print("COMPLETE")
sys.exit(0) sys.exit(0)
......
# intersect.py
#
# Demonstrate how Shapely can be used to analyze and plot the intersection of
# a trajectory and regions in space.
from functools import partial
import random
import pylab
from shapely.geometry import LineString, Point
from shapely.ops import cascaded_union
# Build patches as in dissolved.py
r = partial(random.uniform, -20.0, 20.0)
points = [Point(r(), r()) for i in range(100)]
spots = [p.buffer(2.5) for p in points]
patches = cascaded_union(spots)
# Represent the following geolocation parameters
#
# initial position: -25, -25
# heading: 45.0
# speed: 50*sqrt(2)
#
# as a line
vector = LineString(((-25.0, -25.0), (25.0, 25.0)))
# Find intercepted and missed patches. List the former so we can count them
# later
intercepts = [patch for patch in patches.geoms if vector.intersects(patch)]
misses = (patch for patch in patches.geoms if not vector.intersects(patch))
# Plot the intersection
intersection = vector.intersection(patches)
print(intersection)
assert intersection.geom_type in ['MultiLineString']
if __name__ == "__main__":
# Illustrate the results using matplotlib's pylab interface
pylab.figure(num=None, figsize=(4, 4), dpi=180)
# Plot the misses
for spot in misses:
x, y = spot.exterior.xy
pylab.fill(x, y, color='#cccccc', aa=True)
pylab.plot(x, y, color='#999999', aa=True, lw=1.0)
# Do the same for the holes of the patch
for hole in spot.interiors:
x, y = hole.xy
pylab.fill(x, y, color='#ffffff', aa=True)
pylab.plot(x, y, color='#999999', aa=True, lw=1.0)
# Plot the intercepts
for spot in intercepts:
x, y = spot.exterior.xy
pylab.fill(x, y, color='red', alpha=0.25, aa=True)
pylab.plot(x, y, color='red', alpha=0.5, aa=True, lw=1.0)
# Do the same for the holes of the patch
for hole in spot.interiors:
x, y = hole.xy
pylab.fill(x, y, color='#ffffff', aa=True)
pylab.plot(x, y, color='red', alpha=0.5, aa=True, lw=1.0)
# Draw the projected trajectory
pylab.arrow(-25, -25, 50, 50, color='#999999', aa=True,
head_width=1.0, head_length=1.0)
for segment in intersection.geoms:
x, y = segment.xy
pylab.plot(x, y, color='red', aa=True, lw=1.5)
# Write the number of patches and the total patch area to the figure
pylab.text(-28, 25,
"Patches: %d/%d (%d), total length: %.1f" \
% (len(intercepts), len(patches.geoms),
len(intersection.geoms), intersection.length))
pylab.savefig('intersect.png')
\ No newline at end of file
#!/usr/bin/python
# coding: utf-8
import sys
import pygeoj
import json
from traceback import print_exc
from shapely.geometry import LineString
from shapely.strtree import STRtree
def main(method):
try:
with open('coastline_HELCOM_82214.dat',"r") as coastlinefile:
lines = coastlinefile.readlines()
section = []
sections = []
for line in lines:
if '-999999' in line:
if len(section) > 0:
sections.append(LineString(section))
section = []
continue
split = line.split(' ')
section.append([float(split[2]),float(split[1])])
tree = STRtree(sections)
queryline = LineString([[57.5,17.5],[57.5,18.5]])
matches = tree.query(queryline);
for match in matches:
if match.intersects(queryline):
print(match.intersection(queryline))
# coastfile = ogr.Open("stwcoast.shp")
# coastshape = coastfile.GetLayer(0)
# #first feature of the shapefile
# feature = coastshape.GetFeature(0)
# first = feature.ExportToJson()
# print first # (GeoJSON format)
# lines = gp.GeoDataFrame.from_file('stwcoast.shp')
# intersections= gp.sjoin(poly, lines, how="inner", op='intersects')
# points = [LineString([(i, i), (i+1, i+1)]) for i in range(10)]
# tree = STRtree(points)
# query = tree.query(LineString([(0, 0), (0.5, 0.5)]))
# print([q.wkt for q in query])
sys.exit(0)
except Exception as error:
print_exc()
sys.exit(1)
'''
To run in terminal call with python ex.py
'''
if __name__ == "__main__":
method = ''
if (len(sys.argv) > 1):
method = sys.argv[1]
if (len(sys.argv) > 2):
selected = int(sys.argv[2])
main(method)
\ No newline at end of file
...@@ -19,10 +19,15 @@ ...@@ -19,10 +19,15 @@
import math import math
import json import json
import random import random
import time
from output import Output from output import Output
from output import Particletrack from output import Particletrack
from output import Cloudtrack from output import Cloudtrack
from shapely import geometry from shapely import geometry
from shapely.geometry import LineString
from shapely.strtree import STRtree
from shapely.geometry import MultiLineString
class Model(object): class Model(object):
...@@ -86,7 +91,7 @@ class Model(object): ...@@ -86,7 +91,7 @@ class Model(object):
def createOutlet(self,geom): def createOutlet(self,geom):
nrOfParticles = 300 nrOfParticles = 100
center = geom.centroid center = geom.centroid
points = [] points = []
depth = [] depth = []
...@@ -135,10 +140,68 @@ class Model(object): ...@@ -135,10 +140,68 @@ class Model(object):
return featurecollection return featurecollection
def onlandPoints(self, points, strtree): def onlandPoints(self, beforepoints, nowpoints, matches):
result = list(points) active = []
print("On land calculated") deactive = []
return result 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): def calculateRadius(self, pDist, points):
radius = [pDist] * len(points) radius = [pDist] * len(points)
...@@ -214,22 +277,31 @@ class Model(object): ...@@ -214,22 +277,31 @@ class Model(object):
featurecollection = self.createFeatureCollection() featurecollection = self.createFeatureCollection()
features = featurecollection['features'] features = featurecollection['features']
time = exercisefeature['properties']['time'] timearray = exercisefeature['properties']['time']
properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms)) properties = self.createProperties(1, starttime, multipoint.centroid, [0] * len(multipoint.geoms), [2] * len(multipoint.geoms))
features.append(self.createFeature(multipoint, properties)) features.append(self.createFeature(multipoint, properties))
mp = multipoint mp = multipoint
deactivepoints = []
linestring = geometry.shape(exercisefeature['geometry']) linestring = geometry.shape(exercisefeature['geometry'])
latlngpoints = [[latlng[1], latlng[0]] for latlng in linestring.coords]
matches = strtree.query(LineString(latlngpoints))
print('matches: '+str(matches))
for i, coord in enumerate(linestring.coords): for i, coord in enumerate(linestring.coords):
points = self.centerPoints(coord, mp.centroid, mp.geoms) points = self.centerPoints(coord, mp.centroid, mp.geoms)
displacedpoints = self.displacePoints(pDist, points) displacedpoints = self.displacePoints(pDist, points)
# onlandpoints = self.onlandPoints(displacedpoints,strtree) start = time.time()
mp = geometry.MultiPoint(displacedpoints) activepoints = []
activepoints, onlandpoints = self.onlandPoints(mp.geoms, displacedpoints, matches)
end = time.time()
print('onland: '+str(end-start))
deactivepoints += onlandpoints
mp = geometry.MultiPoint(activepoints+deactivepoints)
level = [0] * len(mp.geoms) level = [0] * len(mp.geoms)
category = [2] * len(mp.geoms) category = [2] * len(mp.geoms)
properties = self.createProperties((i+2), time[i], mp.centroid, level, category) properties = self.createProperties((i+2), timearray[i], mp.centroid, level, category)
features.append(self.createFeature(mp, properties)) features.append(self.createFeature(mp, properties))
if (i % 10 == 0): # if (i % 10 == 0):
Cloudtrack.write(featurecollection) # Cloudtrack.write(featurecollection)
Particletrack.write(featurecollection) # Particletrack.write(featurecollection)
return featurecollection return featurecollection
{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"simulation":{"forwardCalculation":true,"startDate":1513080000000,"stopDate":1513123200000,"oilclass":"Oil classes","substance":"Medium oils (100-1000 cSt)","fresh":true,"amount":"10","amountUnit":"m3","uncertainty":false,"mode":"Normal"},"primary":"PADM","type":"exercise","subtype":"instant","meanLat":57.532041723896384,"meanLng":19.761657714843754,"observation":[{"id":null,"provider":"USER","providerDataRef":null,"providerImgRef":null,"providerType":null,"type":"MANUAL","date":1513080000000}],"domain":"HELCOM","model":"NEMO"},"geometry":{"type":"Polygon","coordinates":[[[19.720458984375004,57.53057086019259],[19.731445312500004,57.54826058225394],[19.775390625,57.54531289147553],[19.802856445312504,57.53057086019259],[19.775390625,57.51582286553883],[19.736938476562504,57.51582286553883],[19.720458984375004,57.51582286553883],[19.720458984375004,57.53057086019259]]]}},{"type":"Feature","properties":{"auxiliary":"exercise","time":[1513080900000,1513081800000,1513082700000,1513083600000,1513084500000,1513085400000,1513086300000,1513087200000,1513088100000,1513089000000,1513089900000,1513090800000,1513091700000,1513092600000,1513093500000,1513094400000,1513095300000,1513096200000,1513097100000,1513098000000,1513098900000,1513099800000,1513100700000,1513101600000,1513102500000,1513103400000,1513104300000,1513105200000,1513106100000,1513107000000,1513107900000,1513108800000,1513109700000,1513110600000,1513111500000,1513112400000,1513113300000,1513114200000,1513115100000,1513116000000,1513116900000,1513117800000,1513118700000,1513119600000,1513120500000,1513121400000,1513122300000,1513123200000]},"geometry":{"type":"LineString","coordinates":[[19.7698974609375,57.53057086019259],[19.778137207031254,57.52614708801571],[19.79049682617188,57.520247890174275],[19.800109863281254,57.51656040692372],[19.811096191406254,57.51139730407089],[19.81796264648438,57.50844663149949],[19.827575683593754,57.50475795522328],[19.83993530273438,57.50033105150468],[19.85641479492188,57.49590361081804],[19.868774414062504,57.492951685358975],[19.882507324218754,57.492213666700756],[19.892120361328125,57.492213666700756],[19.90310668945313,57.492951685358975],[19.91134643554688,57.49516565182901],[19.91958618164063,57.49885529760746],[19.92919921875,57.50254457048255],[19.93881225585938,57.507708926072574],[19.94705200195313,57.512872550877304],[19.953918457031254,57.51729793339678],[19.960784912109375,57.52098534209085],[19.967651367187504,57.52540974047066],[19.973144531250004,57.5290963291026],[19.98001098632813,57.53204533164888],[19.99099731445313,57.534256927023264],[20.00473022460938,57.535731249401564],[20.014343261718754,57.53720551215121],[20.02395629882813,57.537942621165705],[20.03768920898438,57.537942621165705],[20.05142211914063,57.537942621165705],[20.061035156250004,57.537942621165705],[20.0775146484375,57.53646838822989],[20.08712768554688,57.534256927023264],[20.09811401367188,57.52983360210189],[20.10635375976563,57.52540974047066],[20.11184692382813,57.52098534209085],[20.115966796875004,57.51582286553883],[20.118713378906254,57.51139730407089],[20.120086669921875,57.5069712057317],[20.120086669921875,57.498117398284776],[20.120086669921875,57.492213666700756],[20.12283325195313,57.48704711838609],[20.12283325195313,57.48187983904036],[20.11734008789063,57.47671182860216],[20.109100341796875,57.47301994493327],[20.102233886718754,57.469327688204295],[20.096740722656254,57.46489648765881],[20.091247558593754,57.45972607462754],[20.08712768554688,57.45307732607687]]}}]} {"type":"FeatureCollection","features":[{"type":"Feature","properties":{"simulation":{"forwardCalculation":true,"startDate":1513152000000,"stopDate":1513195200000,"oilclass":"Oil classes","substance":"Light oils (0-100 cSt)","fresh":true,"amount":"10","amountUnit":"m3","uncertainty":false,"mode":"Normal"},"primary":"PADM","type":"exercise","subtype":"instant","meanLat":60.00310586337268,"meanLng":27.01194763183594,"observation":[{"id":null,"provider":"USER","providerDataRef":null,"providerImgRef":null,"providerType":null,"type":"MANUAL","date":1513152000000}],"domain":"HELCOM","model":"NEMO"},"geometry":{"type":"Point","coordinates":[27.01194763183594,60.00310586337268]}},{"type":"Feature","properties":{"auxiliary":"exercise","time":[1513152900000,1513153800000,1513154700000,1513155600000,1513156500000,1513157400000,1513158300000,1513159200000,1513160100000,1513161000000,1513161900000,1513162800000,1513163700000,1513164600000,1513165500000,1513166400000,1513167300000,1513168200000,1513169100000,1513170000000,1513170900000,1513171800000,1513172700000,1513173600000,1513174500000,1513175400000,1513176300000,1513177200000,1513178100000,1513179000000,1513179900000,1513180800000,1513181700000,1513182600000,1513183500000,1513184400000,1513185300000,1513186200000,1513187100000,1513188000000,1513188900000,1513189800000,1513190700000,1513191600000,1513192500000,1513193400000,1513194300000,1513195200000]},"geometry":{"type":"LineString","coordinates":[[27.01126098632813,60.00336333025966],[27.01023101806641,60.00362079514246],[27.009372711181644,60.003964078535255],[27.008342742919925,60.004307358365125],[27.007484436035156,60.00465063463208],[27.00645446777344,60.005165542352074],[27.00593948364258,60.00568044205568],[27.00525283813477,60.0061095190184],[27.00473785400391,60.006624404025445],[27.004051208496097,60.007225093068755],[27.003364562988285,60.00765414999103],[27.002849578857425,60.00816901094978],[27.002506256103516,60.00876967193704],[27.00216293334961,60.009284515528016],[27.001819610595703,60.009885156253155],[27.001647949218754,60.01057158943673],[27.0014762878418,60.01117220678301],[27.001304626464847,60.01168701296785],[27.001132965087894,60.01228761005297],[27.001132965087894,60.012973993362905],[27.001132965087894,60.01357456707032],[27.001132965087894,60.014260923663215],[27.001132965087894,60.01486147399353],[27.001132965087894,60.01546201341474],[27.001132965087894,60.016062541926956],[27.001132965087894,60.01683463398483],[27.001132965087894,60.01760670801],[27.001304626464847,60.01837876400275],[27.001304626464847,60.019150801963356],[27.0014762878418,60.019751263466226],[27.001647949218754,60.02060904669113],[27.001819610595703,60.02112370594082],[27.001991271972656,60.021809905807785],[27.001991271972656,60.02241031900469],[27.001991271972656,60.023096492159375],[27.001991271972656,60.02395418856897],[27.001647949218754,60.02472609630581],[27.0014762878418,60.02524069144713],[27.00096130371094,60.02609833220896],[27.000789642333988,60.026612905981985],[27.00044631958008,60.02738475161725],[27.00044631958008,60.02807082149105],[27.00044631958008,60.0286711209454],[27.00044631958008,60.029357164110294],[27.00044631958008,60.030043193030735],[27.00027465820313,60.03055770537328],[26.999931335449222,60.031329458864015],[26.99975967407227,60.031929699115615]]}}]}
\ No newline at end of file \ No newline at end of file
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment