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.

marked as duplicate by user2856, MappaGnosis, Mapperz Nov 12 '15 at 18:45

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • Can you provide a small subset for testing as well as the desired outcome? This would make the explanation using Python code a lot easier. – Kersten Nov 11 '15 at 7:48
  • Polygons: [link] (drive.google.com/…) Points: [link] (drive.google.com/…) Desired output: shapefile containing points that fall within polygons. Or simply a count of the number of points falling within polygons. – geoeye Nov 12 '15 at 22:33
  • And where does the raster mentioned in the question come into play? Your original question sounds like zonal statistics/extract raster values below points. – Kersten Nov 13 '15 at 8:31

Browse other questions tagged or ask your own question.