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 <[email protected]> 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 > [email protected] > https://lists.osgeo.org/mailman/listinfo/gdal-dev > -- Sean Gillies
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
