[gdal-dev] Querying RasterBands near edges

Andrew C Aitchison andrew at aitchison.me.uk
Mon Dec 7 13:45:44 PST 2020


On Mon, 7 Dec 2020, Felix wrote:

> 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.

Unless your dataset covers the whole world, I don't understand why would 
wld want such a wrap-around for a geographical application.
Even then you would wrap at the left/right border but not the poles 
(top/bottom).

If you are processing whole world data, could you generate a second raster 
with the left and right halves swapped and split the points into two sets ?
Then provided you use the correct image for each set your radius wont 
cross the edge of the raster.

-- 
Andrew C. Aitchison					Kendal, UK
 			andrew at aitchison.me.uk


More information about the gdal-dev mailing list