[gdal-dev] Overviews are not taken into account while reading with specified resampling method
Denis Rykov
rykovd at gmail.com
Thu Aug 27 03:41:23 PDT 2020
I was able to reproduce this issue with pure GDAL. When you read data with
boundless=True in rasterio it creates an intermediate VRT file. This is the
example of file that being created in my case:
<?xml version="1.0" encoding="UTF-8"?>
<VRTDataset rasterXSize="40961" rasterYSize="139265">
<SRS>PROJCS["WGS 84 / UTM zone 47N",GEOGCS["WGS
84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",99],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32647"]]</SRS>
<GeoTransform>443070.0,1.0,0.0,4366312.0,0.0,-1.0</GeoTransform>
<VRTRasterBand band="1" dataType="Byte">
<NoDataValue>0.0</NoDataValue>
<ColorInterp>Red</ColorInterp>
<ComplexSource>
<SourceFilename relativeToVRT="1" shared="0">dummy.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="1" yOff="0" ySize="1" />
<DstRect xOff="0" xSize="1" yOff="0" ySize="1" />
<ScaleRatio>0</ScaleRatio>
<ScaleOffset>0.0</ScaleOffset>
</ComplexSource>
<ComplexSource>
<SourceFilename relativeToVRT="0"
shared="0">/vsicurl/https://*.vrt</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="40960" yOff="0" ySize="139264" />
<DstRect xOff="-6961.0" xSize="40960.0" yOff="-105176.0" ySize="139264.0" />
<NODATA>0.0</NODATA>
<OpenOptions /> </ComplexSource>
</VRTRasterBand>
<VRTRasterBand band="2" dataType="Byte">
<NoDataValue>0.0</NoDataValue>
<ColorInterp>Green</ColorInterp>
<ComplexSource>
<SourceFilename relativeToVRT="1" shared="0">dummy.tif</SourceFilename>
<SourceBand>2</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="1" yOff="0" ySize="1" />
<DstRect xOff="0" xSize="1" yOff="0" ySize="1" />
<ScaleRatio>0</ScaleRatio>
<ScaleOffset>0.0</ScaleOffset>
</ComplexSource>
<ComplexSource>
<SourceFilename relativeToVRT="0"
shared="0">/vsicurl/https://*.vrt</SourceFilename>
<SourceBand>2</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="40960" yOff="0" ySize="139264" />
<DstRect xOff="-6961.0" xSize="40960.0" yOff="-105176.0" ySize="139264.0" />
<NODATA>0.0</NODATA>
<OpenOptions /> </ComplexSource>
</VRTRasterBand>
<VRTRasterBand band="3" dataType="Byte">
<NoDataValue>0.0</NoDataValue>
<ColorInterp>Blue</ColorInterp>
<ComplexSource>
<SourceFilename relativeToVRT="1" shared="0">dummy.tif</SourceFilename>
<SourceBand>3</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="1" yOff="0" ySize="1" />
<DstRect xOff="0" xSize="1" yOff="0" ySize="1" />
<ScaleRatio>0</ScaleRatio>
<ScaleOffset>0.0</ScaleOffset>
</ComplexSource>
<ComplexSource>
<SourceFilename relativeToVRT="0"
shared="0">/vsicurl/https://*.vrt</SourceFilename>
<SourceBand>3</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="40960" yOff="0" ySize="139264" />
<DstRect xOff="-6961.0" xSize="40960.0" yOff="-105176.0" ySize="139264.0" />
<NODATA>0.0</NODATA>
<OpenOptions /> </ComplexSource>
</VRTRasterBand>
<MaskBand>
<VRTRasterBand dataType="Byte">
<SimpleSource>
<SourceFilename relativeToVRT="0"
shared="0">/vsicurl/https://*.vrt</SourceFilename>
<SourceBand>mask,1</SourceBand>
<SourceProperties BlockXSize="128" BlockYSize="128"
RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
<SrcRect xOff="0" xSize="40960" yOff="0" ySize="139264" />
<DstRect xOff="-6961.0" xSize="40960" yOff="-105176.0"
ySize="139264" /> </SimpleSource>
</VRTRasterBand>
</MaskBand>
</VRTDataset>
If I read data from this file with gdal using
resample_alg=gdalconst.GRIORA_NearestNeighbour then GDAL takes into
account *.vrt.ovr file and sends very few HTTP requests to the server
(~30):
ds = gdal.OpenEx("/tmp/rasterio.vrt")
image = ds.ReadAsArray(xoff=0, yoff=0, xsize=5671, ysize=5648,
buf_xsize=383, buf_ysize=385,
resample_alg=gdalconst.GRIORA_NearestNeighbour
If I do the same but using resample_alg=gdalconst.GRIORA_Cubic then GDAL
sends a huge amount of requests to the server (~1k) because overviews are
not used:
ds = gdal.OpenEx("/tmp/rasterio.vrt")
image = ds.ReadAsArray(xoff=0, yoff=0, xsize=5671, ysize=5648,
buf_xsize=383, buf_ysize=385, resample_alg=gdalconst.GRIORA_Cubic
Is it expected or might there be something wrong with that VRT file? Thanks
in advance for any help.
On Wed, Aug 26, 2020 at 6:18 PM Denis Rykov <rykovd at gmail.com> wrote:
> I have remote *.vrt raster and *.vrt.ovr accessible through HTTP. When I
> run the following script with rasterio:
>
> with rasterio.open("http://*.vrt"") as src:
> image = src.read(indexes=[1, 2, 3], **{
> "window": Window(col_off=6961, row_off=105176, width=5671, height=5648),
> "resampling": Resampling.cubic,
> "boundless": True,
> "out_shape": (3, 383, 385),
> "masked": True
> })
>
> depending on "resampling" algorithm GDAL sends different amounts of
> requests to the server. In the case of "cubic" it doesn't take into account
> overviews and sends requests directly to *.tif files (900 in my case). In
> case of "nearest" everything is ok (only 60 requests, *.vrt.ovr is taken
> into account).
>
> Does GDAL check the resampling algorithm of overviews and in case it
> differs from the option specified in read() method they are bypassed or it
> works differently?
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200827/7b0fcdcd/attachment-0001.html>
More information about the gdal-dev
mailing list