[gdal-dev] Querying RasterBands near edges

[Sean Gillies]

Hi Felix,

The rasterio Python API has a concept of a "boundless" read (see boundless
keyword under
https://rasterio.readthedocs.io/en/latest/api/rasterio._io.html#rasterio._io.DatasetReaderBase.read)
which sounds similar to what you're looking for. This is implemented using
a VRT and you might be able to do the same kind of thing with the GDAL
Python bindings.

1. Create a VRT with ample extra space around your raster and with a
suitable background value
2. Read from the VRT instead of the raster

You could even tile the VRT with multiple virtual copies of your raster if
you wanted a true wrap-around.

Hope this helps!

On Mon, Dec 7, 2020 at 1:59 PM Felix <felix.divo at sailingteam.tu-darmstadt.de>
wrote:

> 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.
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>


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