I have several polygons which are actually a union of points. I want a reasonably quick way of removing the points which are not land (but on river, lake, ocean, etc).
So I far I have come up with the below method which takes me from the left picture to the right one:
import fiona
from shapely.geometry import Point, shape
# Shape files (Small Lakes Europe, Lakes, Coast)
le_loc = ".../water-polygons-split-4326/ne_10m_lakes_europe.shp"
lw_loc = ".../water-polygons-split-4326/ne_10m_lakes.shp"
c_loc = ".../water-polygons-split-4326/water_polygons.shp"
print("Opening Shape-Files")
le_poly = [shape(pol['geometry']) for pol in fiona.open(le_loc)]
lw_poly = [shape(pol['geometry']) for pol in fiona.open(lw_loc)]
c_poly = [shape(pol['geometry']) for pol in fiona.open(c_loc)]
def point_in_water(lat, lon):
xy_point = Point(lon, lat)
for one in [c_poly, lw_poly, le_poly]:
for shape in one:
if shape.contains(xy_point):
return True
return False
# Test (true)
print(point_in_water(46.268408, 6.180437))
Then in a loop I test my points:
with open(in_csv) as f:
for x in csv.reader(f):
# Lat, Lng, Mins, Mode
if not point_in_water(float(x[0]), float(x[1])):
coords.append([x[0], x[1])
I use three shape-files which are OK for my purpose (the lakes one is a bit rough): coast, lakes, small lakes.
However, for 10,000 points the code is a bit slow (and I have around 30 files so 300,000 points to examine).
I was wondering if any of the following would be possible:
1) I am looping through shapes and checking shape.contains(point) instead of looping through points and checking point.within(shape) -> I wasn't sure if there would be an improvement?
2) Using a spatial-index would probably speed this up, however I don't think RTree works anymore with Python 3.4
3) Perhaps there is a faster function (a rough contains) which just checks bounds and I can use that as step one and then contains as step 2.
4) Instead of looping through points, is there a way to vectorise and pass all at once?
5) Perhaps it would be faster to convert the shapely Polygons to a path so that I can use matpotlib's path.contains_point?
6) Finally, I realise that I should probably use a Mercator projection for the point in poly test, but a rough cut is fine for me (it's not like the water shapefile is super accurate anyway).
Thank you