Skip to content

Instantly share code, notes, and snippets.

@richpsharp
Created May 13, 2020 04:41
Show Gist options
  • Save richpsharp/cac0d490c0b9ee5ce0f64015a0d14b4b to your computer and use it in GitHub Desktop.
Save richpsharp/cac0d490c0b9ee5ce0f64015a0d14b4b to your computer and use it in GitHub Desktop.
Example script for testing points contained in a set of polygons
"""Example of determing if points are contained in polygon."""
from osgeo import gdal
import shapely.geometry
point_vector_path = 'dams.shp'
polygon_vector_path = 'watersheds_hotels_resorts_DdS.shp'
polygon_id_field = 'Name'
point_id_field = 'desc_short'
if __name__ == "__main__":
polygon_vector = gdal.OpenEx(polygon_vector_path, gdal.OF_VECTOR)
polygon_layer = polygon_vector.GetLayer()
point_vector = gdal.OpenEx(point_vector_path, gdal.OF_VECTOR)
point_layer = point_vector.GetLayer()
for polygon_feature in polygon_layer:
polygon_geom = polygon_feature.GetGeometryRef()
point_layer.ResetReading()
print(
f"testing polygon: {polygon_feature.GetField(polygon_id_field)}")
for point in point_layer:
point_geom = point.GetGeometryRef()
if polygon_geom.Intersects(point_geom):
print(
f"\t{point.GetField(point_id_field)} intersects "
f"{polygon_feature.GetField(polygon_id_field)}")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment