This question already has an answer here:
I have a raster with one band and I would like to query a set of overlying points (.shp) to determine how many points as well as which points points fall under a range of certain pixel values associated with the raster.
This is straightforward in ArcGIS but I want to use GDAL/OGR/Python to automate this process.
I can think of a few ways to potentially accomplish this:
(1) Join the raster values to the points, add a new field with the raster values and then query the resulting new_points.shp file.
(2) Write an effective SQL statement that queries the points and the raster to write out a new file with only those points that correspond to the value range I'm interested in. Potentially using a VRT file
(3) Convert the raster to a shapefile using gdal_polygonize (very time intensive!) and then use an SQL statement along the lines of
SELECT COUNT(*) FROM points, new_poly WHERE ST_INTERSECTS(new_poly.geometry,points.geometry) -dialect SQLITE input.vrt again using a vrt table.
Number (2) seems like the best and most efficient option, but I am unfamiliar with querying both raster and vector data simultaneously. Is there a standard or accepted way to do this?
EDIT: Another option seems to be the
gdallocationinfo utility, but it seems that you can only pass lat/long coordinates and not shapefiles.