[gdal-dev] Overviews are not taken into account while reading with specified resampling method

Denis Rykov rykovd at gmail.com
Thu Aug 27 06:08:02 PDT 2020


I found the culprit. If remove this section from each band definition in
VRT file then everything works fine:

<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

Rasterio's code that builds this code was altered 2 years ago after
conversation in OSGeo/gdal#1135 <https://github.com/OSGeo/gdal/issues/1135>.
But it has this undesirable consequence I'm experiencing. Any ideas how to
overcome that?

On Thu, Aug 27, 2020 at 12:41 PM Denis Rykov <rykovd at gmail.com> wrote:

> 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/dac058ec/attachment-0001.html>


More information about the gdal-dev mailing list