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
