2

I have a EPSG:4326 raster file with pixel size of 15 arcseconds, and the longitude and latitude of a point.

dataset = rasterio.open(url)
lat = 5.535456  
lon = 7.531536

I know I can retrieve the pixel value of the point with the code below.

row, col = dataset.index(lon, lat)
band = dataset.read(1)
band[row][col]

But how can I retrieve the pixel value of a point X km away to the north of the current point? My ultimate intention is to aggregate over a circle of X km radius, whose center is a given point like the one above.

1 Answer 1

6

If you're ultimately working toward a polygon (circle) intersection with your raster, there are built-in methods for calculating a numpy mask, given a polygon.

You can use this question to see how Shapely can be used to get a circle of defined radius from a centre point (you will need to convert to a projected coordinate system).

Then, use rasterio.features.geometry_mask to calcualate a numpy mask given a polygon geometry. (Shapely's abstract geometry interface is widely accepted in an ecosystem of spatial Python libraries, which includes Rasterio.)

Once you have the mask for the ndarray, you can calculate various statistics on the masked array using numpy methods, e.g. numpy.ma.mean.

1
  • Is there any way that I can avoid using numpy masked arrays? I want to perform this operation over list of longitudes and latitudes, and numpy.ma.mean is too slow. Commented Jan 18, 2022 at 15:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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