[gdal-dev] Querying RasterBands near edges

[Felix]

Hey,

Our application requires us to fetch for many pairs of (latitude, 
longitude) all existing values within some radius from an existing 
raster dataset. Is it possible to read values from a GDALRasterBand 
given some general bounding box, even if it is at the edge of a dataset?

One can easily convert between bounding boxes in geographical 
coordinates and pixels using Dataset::GetGeoTransform() and 
InvGeoTransform(), but GDALRasterBand/GDALMDArray/the Python bindings 
seems to not have a method for being queried by it. Of course, simple 
boxes can be read e.g. using the Python methods 
Band::ReadAsArray/Band::ReadRaster. But this does not work once queries 
get to the boundaries, as all of these methods do not support "wrapping 
around" to the other side of the dataset. As an example:

    Dataset: 1000 x 1000 pixels
    Bounding box in pixels, given as two corner points in (x, y) pixels:
    (950, 950) and (50, 50)
    This should query all 50 by 50 pixel corners of the dataset and
    result in a 100 by 100 pixel result
    It arises when a point at the edge of a dataset is queried with a
    radius of 50 pixels
    However, numpy and the above methods do not allow to express such
    wrap-around indexing

Is there some method that I'm overlooking or is this intentionally 
something that users should implement for themselfs? It feels like such 
indexing would be required in a lot of other places too.

I've posted here last week about some suggestions on how to query raster 
bands at non-structured points. Thanks for the useful replies, Christoph 
& Joaquim! However, it seems like we have to code that ourselfs, thus 
the above digging.

Cheers
Felix D.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20201207/3c98e13c/attachment.html>