Find points that intersect a user-defined polygon in GDAL

Question

I'm trying to do something that I don't know is possible in GDAL/Python.

I have a list, L, of pre-defined lat/long points that make up a uniform grid.

A user will send me a set of points that make up the vertices of a polygon. This polygon will usually be irregular, not square. Any number of points could make up the polygon vertices, but it will likely not exceed 10 or so.

I would like to: 1) Create a polygon, P, based on the user-defined vertices. 2) Find the lat/long points in L that intersect polygon P.

Is this possible? I don't normally work with vectors in GDAL, so I'm not familiar with the best way to do this.

Answer 1

You don't need Gdal/Python here. It is easier to use Shapely

1) Transform the list of points to a shapely geometry

list = [(0.0, 0.0), (0.0, 1.0), (0.0, 2.0), (0.0, 3.0), (0.0, 4.0), (0.0, 5.0), (0.0, 6.0), (0.0, 7.0), (0.0, 8.0), (0.0, 9.0), (1.0, 0.0), (1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0), (1.0, 5.0), (1.0, 6.0), (1.0, 7.0), (1.0, 8.0), (1.0, 9.0), (2.0, 0.0), (2.0, 1.0), (2.0, 2.0), (2.0, 3.0), (2.0, 4.0), (2.0, 5.0), (2.0, 6.0), (2.0, 7.0), (2.0, 8.0), (2.0, 9.0), (3.0, 0.0), (3.0, 1.0), (3.0, 2.0), (3.0, 3.0), (3.0, 4.0), (3.0, 5.0), (3.0, 6.0), (3.0, 7.0), (3.0, 8.0), (3.0, 9.0), (4.0, 0.0), (4.0, 1.0), (4.0, 2.0), (4.0, 3.0), (4.0, 4.0), (4.0, 5.0), (4.0, 6.0), (4.0, 7.0), (4.0, 8.0), (4.0, 9.0), (5.0, 0.0), (5.0, 1.0), (5.0, 2.0), (5.0, 3.0), (5.0, 4.0), (5.0, 5.0), (5.0, 6.0), (5.0, 7.0), (5.0, 8.0), (5.0, 9.0), (6.0, 0.0), (6.0, 1.0), (6.0, 2.0), (6.0, 3.0), (6.0, 4.0), (6.0, 5.0), (6.0, 6.0), (6.0, 7.0), (6.0, 8.0), (6.0, 9.0), (7.0, 0.0), (7.0, 1.0), (7.0, 2.0), (7.0, 3.0), (7.0, 4.0), (7.0, 5.0), (7.0, 6.0), (7.0, 7.0), (7.0, 8.0), (7.0, 9.0), (8.0, 0.0), (8.0, 1.0), (8.0, 2.0), (8.0, 3.0), (8.0, 4.0), (8.0, 5.0), (8.0, 6.0), (8.0, 7.0), (8.0, 8.0), (8.0, 9.0), (9.0, 0.0), (9.0, 1.0), (9.0, 2.0), (9.0, 3.0), (9.0, 4.0), (9.0, 5.0), (9.0, 6.0), (9.0, 7.0), (9.0, 8.0), (9.0, 9.0)]
from shapely.geometry import Point, Polygon, MultiPoint
grid = MultiPoint([Point(i,j) for i,j in list])

2) A set of points that make up the vertices of a polygon

pts = [(0.89861841430257, 5.608289514184579), (2.2896616734212287, 6.420009426460781), (3.028966167342123, 5.420676653157543), (3.9131014979736314, 5.261362732775868), (3.5799905735392192, 3.391043259118658), (3.072415418355307, 3.434492510131842), (1.5944736572102807, 3.666889075565587), (0.00032006681885698995, 5.434839670009748), (0.89861841430257, 5.608289514184579)]

3) Create a polygon from this list of points

polygon = Polygon(pts)

4) Use the Shapely predicates

pts =  [pt for pt in points if pt.within(polygon )]
pts2 = [pt for pt in points if pt.intersects(polygon )]
pts3 = [pt for pt in points if pt.touches(polygon )]
for pt in pts:
   print pt,
POINT (1 5) POINT (2 4) POINT (2 5) POINT (2 6) POINT (3 4) POINT (3 5)
for pt in pts2:
   print pt,
POINT (1 5) POINT (2 4) POINT (2 5) POINT (2 6) POINT (3 4) POINT (3 5)
for pt in pts3:
   print pt,
Nothing

5) If you want, save the result with Fiona

import fiona
from shapely.geometry import mapping
# schema of the shapefile
schema = {'geometry': 'Point', 'properties': {'ido':'int:10'}}
with fiona.open('result.shp', 'w', 'ESRI Shapefile', schema) as ouput:
    for index, pts in enumerate(pts):
      output.write({'geometry': mapping(pts),'properties': {'id':index}})

Answer 2

For your first question, I believe this previous asked question will help you: Why does this simple Python OGR code create an empty polygon?

Here is a python GDAL cookbook (a few years old, but still applicable) GDAL/OGR python Cookbook: https://pcjericks.github.io/py-gdalogr-cookbook/

For your second question, I believe that you want to set a spatial filter. Reference: https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#spatial-filter

There is simpler python library called fiona that may also be up your alley.
http://toblerity.org/fiona/manual.html

This isn't a full answer, but I at least wanted to get you in the right direction.