# Find points that intersect a user-defined polygon in GDAL

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.

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}})
``````
• This is a great solution as long as the grid you are intersecting with is not too big. I'm trying to intersect a world 1km grid with a circle shapefile, and any operations immediately kill my python session. Jan 9 at 9:05

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.