<div dir="ltr">Hi Felix,<br><div><br></div><div>The rasterio Python API has a concept of a "boundless" read (see boundless keyword under <a href="https://rasterio.readthedocs.io/en/latest/api/rasterio._io.html#rasterio._io.DatasetReaderBase.read">https://rasterio.readthedocs.io/en/latest/api/rasterio._io.html#rasterio._io.DatasetReaderBase.read</a>) 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.</div><div><br></div><div>1. Create a VRT with ample extra space around your raster and with a suitable background value</div><div>2. Read from the VRT instead of the raster</div><div><br></div><div>You could even tile the VRT with multiple virtual copies of your raster if you wanted a true wrap-around.</div><div><br></div><div>Hope this helps!</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Dec 7, 2020 at 1:59 PM Felix <<a href="mailto:felix.divo@sailingteam.tu-darmstadt.de">felix.divo@sailingteam.tu-darmstadt.de</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
  

    
  
  <div>
    <p>Hey,</p>
    <p>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?</p>
    <p>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:</p>
    <blockquote>
      <p>Dataset: 1000 x 1000 pixels<br>
        Bounding box in pixels, given as two corner points in (x, y)
        pixels: (950, 950) and (50, 50)<br>
        This should query all 50 by 50 pixel corners of the dataset and
        result in a 100 by 100 pixel result<br>
        It arises when a point at the edge of a dataset is queried with
        a radius of 50 pixels<br>
        However, numpy and the above methods do not allow to express
        such wrap-around indexing</p>
    </blockquote>
    <p>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.<br>
    </p>
    <p>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.</p>
    <p>Cheers<br>
      Felix D.<br>
    </p>
  </div>

_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr">Sean Gillies</div></div>