3

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.

5

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])

enter image description here

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)]

enter image description here

3) Create a polygon from this list of points

polygon = Polygon(pts)

enter image description here

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

enter image description here

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}})
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.