Mercurial > hg > nsaunier > traffic-intelligence

comparison python/moving.py @ 781:7c38250ddfc7 dev

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
updated to deal with prepared polygons from shapely, and to extract the same positions from a second trajectory in a polygon (for velocities for example)
author Nicolas Saunier <nicolas.saunier@polymtl.ca>
date Thu, 11 Feb 2016 11:56:50 -0500
parents 1b22d81ef5ff
children 3aa6102ccc12
comparison
equal deleted inserted replaced
780:1b22d81ef5ff 781:7c38250ddfc7
10 from scipy.stats import scoreatpercentile 10 from scipy.stats import scoreatpercentile
11 from scipy.spatial.distance import cdist 11 from scipy.spatial.distance import cdist
12 12
13 try: 13 try:
14 from shapely.geometry import Polygon, Point as shapelyPoint 14 from shapely.geometry import Polygon, Point as shapelyPoint
15 from shapely.prepared import prep 15 from shapely.prepared import prep, PreparedGeometry
16 shapelyAvailable = True 16 shapelyAvailable = True
17 except ImportError: 17 except ImportError:
18 print('Shapely library could not be loaded') 18 print('Shapely library could not be loaded')
19 shapelyAvailable = False 19 shapelyAvailable = False
20 20
347 'Returns the middle of the segment [p1, p2]' 347 'Returns the middle of the segment [p1, p2]'
348 return Point(0.5*p1.x+0.5*p2.x, 0.5*p1.y+0.5*p2.y) 348 return Point(0.5*p1.x+0.5*p2.x, 0.5*p1.y+0.5*p2.y)
349 349
350 if shapelyAvailable: 350 if shapelyAvailable:
351 def pointsInPolygon(points, polygon): 351 def pointsInPolygon(points, polygon):
352 '''Optimized tests of a series of points within (Shapely) polygon ''' 352 '''Optimized tests of a series of points within (Shapely) polygon (not prepared)'''
353 prepared_polygon = prep(polygon) 353 if type(polygon) == PreparedGeometry:
354 prepared_polygon = polygon
355 else:
356 prepared_polygon = prep(polygon)
354 return filter(prepared_polygon.contains, points) 357 return filter(prepared_polygon.contains, points)
355 358
356 # Functions for coordinate transformation 359 # Functions for coordinate transformation
357 # From Paul St-Aubin's PVA tools 360 # From Paul St-Aubin's PVA tools
358 def subsec_spline_dist(splines): 361 def subsec_spline_dist(splines):
931 'Returns the positions very step' 934 'Returns the positions very step'
932 return Trajectory([self.positions[0][::step], 935 return Trajectory([self.positions[0][::step],
933 self.positions[1][::step]]) 936 self.positions[1][::step]])
934 937
935 if shapelyAvailable: 938 if shapelyAvailable:
936 def getTrajectoryInPolygon(self, polygon): 939 def getTrajectoryInPolygon(self, polygon, t2 = None):
937 '''Returns the trajectory built with the set of points inside the (shapely) polygon''' 940 '''Returns the trajectory built with the set of points inside the (shapely) polygon
941 The polygon could be a prepared polygon (faster) from prepared.prep
942
943 t2 is another trajectory (could be velocities)
944 which is filtered based on the first (self) trajectory'''
938 traj = Trajectory() 945 traj = Trajectory()
939 points = [p.asShapely() for p in self] 946 inPolygon = []
940 for p in pointsInPolygon(points, polygon): 947 for x, y in zip(self.positions[0], self.positions[1]):
941 traj.addPositionXY(p.x, p.y) 948 inPolygon.append(polygon.contains(shapelyPoint(x, y)))
942 return traj 949 if inPolygon[-1]:
950 traj.addPositionXY(x, y)
951 traj2 = Trajectory()
952 if t2 is not None:
953 for inp, x, y in zip(inPolygon, t2.positions[0], t2.positions[1]):
954 if inp:
955 traj2.addPositionXY(x, y)
956 return traj, traj2
943 957
944 def proportionInPolygon(self, polygon, minProportion = 0.5): 958 def proportionInPolygon(self, polygon, minProportion = 0.5):
945 pointsIn = pointsInPolygon([p.asShapely() for p in self], polygon) 959 inPolygon = [polygon.contains(shapelyPoint(x, y)) for x, y in zip(self.positions[0], self.positions[1])]
946 lengthThreshold = float(self.length())*minProportion 960 lengthThreshold = float(self.length())*minProportion
947 return len(pointsIn) >= lengthThreshold 961 return sum(inPolygon) >= lengthThreshold
948 else: 962 else:
949 def getTrajectoryInPolygon(self, polygon): 963 def getTrajectoryInPolygon(self, polygon, t2 = None):
950 '''Returns the trajectory built with the set of points inside the polygon 964 '''Returns the trajectory built with the set of points inside the polygon
951 (array of Nx2 coordinates of the polygon vertices)''' 965 (array of Nx2 coordinates of the polygon vertices)'''
952 traj = Trajectory() 966 traj = Trajectory()
967 inPolygon = []
953 for p in self: 968 for p in self:
954 if p.inPolygon(polygon): 969 inPolygon.append(p.inPolygon(polygon))
970 if inPolygon[-1]:
955 traj.addPosition(p) 971 traj.addPosition(p)
956 return traj 972 traj2 = Trajectory()
973 if t2 is not None:
974 for inp, x, y in zip(inPolygon, t2.positions[0], t2.positions[1]):
975 if inp:
976 traj2.addPositionXY(p.x, p.y)
977 return traj, traj2
957 978
958 def proportionInPolygon(self, polygon, minProportion = 0.5): 979 def proportionInPolygon(self, polygon, minProportion = 0.5):
959 pointsInPolygon = [p.inPolygon(polygon) for p in self] 980 inPolygon = [p.inPolygon(polygon) for p in self]
960 lengthThreshold = float(self.length())*minProportion 981 lengthThreshold = float(self.length())*minProportion
961 return len(pointsInPolygon) >= lengthThreshold 982 return sum(inPolygon) >= lengthThreshold
962 983
963 @staticmethod 984 @staticmethod
964 def lcss(t1, t2, lcss): 985 def lcss(t1, t2, lcss):
965 return lcss.compute(t1, t2) 986 return lcss.compute(t1, t2)
966 987