To find which feature the coordinate is near, you first need to build an R-Tree index of bounding boxes or envelope of each feature. A popular library for this is libspatialindex.
Secondly, you would then need to know for each of the matched features from your R-tree, which ones match. GDAL/OGR does have some operations based on GEOS to see if the geometry "contains" the point of interest, and you can extract the field info. See the OGRGeometry Class Reference for Contains, Touches, Within, etc. that can perform the geometry relation operator.
Personally I wouldn't do any of the above, because it would take me too long to figure out the coding, and I know there are faster ways. (I would loading all shapefiles in a PostGIS database, then query the locations using SQL.) However, it is understandable if this GIS functionality needs to be embedded in an existing project without adding more complicated dependencies, such as a database server.
It is a pure Python problem and not an GDAL/OGR problem. But if you want to program in Python,you need to understand the difference between an iterator and an iterable (many, many references on the Web, but one of the best is Loop Like A Native by Ned Batchelder). These are fundamental concepts in Python.
"The important operation on an iterable is iter(), which will return an iterator. And the only operation available on an iterator is next(), which either returns the next value, or raises StopIteration, a special exception that means the iteration is finished." (from Ned Batchelder)
An iterable produces an iterator and an iterator produces a stream of values.
So what's the difference between the for loop and pointsLayer.GetNextFeature() ?
the for loops = iterator and pointsLayer= iterable
for feature in iterable:
statements
so with your example:
from osgeo import ogr
source = ogr.Open('yourpointshapefile.shp')
pointsLayer = source.GetLayer()
for feature in pointsLayer:
geom =feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272022.68669955182, 155404.12013724342, 0.0)
(272904.99338241993, 152881.6706538822, 0.0)
.....
(272796.14718989708, 152075.00336062044, 0.0)
But we can create another type of iterator with iter() and next():
source = ogr.Open('yourpointshapefile.shp')
pointsLayer = source.GetLayer()
iterator = iter(pointsLayer)
# first feature in pointsLayer
feature = iterator.next()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272022.68669955182, 155404.12013724342, 0.0)
# second feature in pointsLayer
feature = iterator.next()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272904.99338241993, 152881.6706538822, 0.0)
....
# end raises StopIteration error to signal that iteration is complete
feature = iterator.next()
Traceback (most recent call last):
...
StopIteration
Unlike the for loop that handles the StopIteration error, you need here another control structure (while loop,if,...)
Thus, what is pointsLayer.GetNextFeature() ?
It is an iterator and you can replace iterator = iter(pointsLayer) and iterator.next() by
feature = pointsLayers.GetNextFeature()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
.....
feature = pointsLayers.GetNextFeature()
....
I hope that I have been able to explain why you should not mix the for loops (iterator) with .GetNextFeature() (another iterator).
Best Answer
You should be able to use the ogr2ogr utility to extract the features you want from your shapefile, using a bounding box.
Something like:
and also by specifying a polygon instead of a bounding box like: