Re: [gdal-dev] Overviews are not taken into account while reading with specified resampling method
Hi Denis, Even, On Thu, Aug 27, 2020 at 8:08 AM Even Rouault wrote: > On jeudi 27 août 2020 15:08:02 CEST Denis Rykov wrote: > > > I found the culprit. If remove this section from each band definition in > > > VRT file then everything works fine: > > > > > > > > > dummy.tif > > > 3 > > > > > RasterXSize="40961" RasterYSize="139265" dataType="Byte" /> > > > > > > > > > 0 > > > 0.0 > > > I'm perhaps missing something, but the above snippet added by RasterIO > doesn't make sense to me. It just sets the pixel at (0,0) to 0. I bet it is > completely useless. It would make more sense to have xSize/ySize of DstRect > to cover the whole raster. And as is set on the > VRTRasterBand, this should be used to fill the target raster anyway. > > > > As far as why overviews of the original VRT aren't used is concerned, I'm > not sure why. I'd have expected it to work, but I must be missing something. > > > > Even > That's a rasterio bug: the background "dummy.tif" source should be . I'm working on a patch for this. This VRT background fill is different from nodata. The scale offset could be different from the nodata value. BTW, It takes advantage of GDAL shortcuts predicated on scale ratio. The "dummy.tif" file is never opened and doesn't even exist. -- Sean Gillies ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Overviews are not taken into account while reading with specified resampling method
In GDAL 3.1 (previous example was done in GDAL 2.2.3) I can read data from that VRT using "gdalconst.GRIORA_Cubic" but still overviews are not used: >>> from osgeo import gdal >>> from osgeo.gdal import gdalconst >>> gdal.__version__ '3.1.0' >>> url = >>> "/vsicurl/https://gist.githubusercontent.com/drnextgis/3cbdbace7b5b8b80c3c6169b109bf9db/raw/7b61af41c698930a7f2da851a0684a15ddeb99a9/rasterio-boundless.vrt; >>> ds = gdal.OpenEx(url) >>> image = ds.ReadAsArray(xoff=0, yoff=0, xsize=64, ysize=64, buf_xsize=10, >>> buf_ysize=10, resample_alg=gdalconst.GRIORA_Cubic) >>> image.shape (3, 10, 10) On Thu, Aug 27, 2020 at 7:30 PM Denis Rykov wrote: > Hi Sean. I patched rasterio as you suggested and intermediate file now > looks like this (I'm trying now with a public dataset): > https://gist.github.com/drnextgis/3cbdbace7b5b8b80c3c6169b109bf9db > > But when I read it with GDAL using non-nearest algorithm I'm getting > the following error: > > >>> from osgeo import gdal>>> from osgeo.gdal import gdalconst>>> url = > >>> "/tmp/rasterio-boundless.vrt">>> ds = gdal.OpenEx(url)>>> image = > >>> ds.ReadAsArray(xoff=0, yoff=0, xsize=64, ysize=64, buf_xsize=10, > >>> buf_ysize=10, resample_alg=gdalconst.GRIORA_Cubic) > ERROR 4: /tmp/dummy.tif: No such file or directory > > but with NearestNeighbour it works without error: > > >>> from osgeo import gdal>>> from osgeo.gdal import gdalconst>>> url = > >>> "/tmp/rasterio-boundless.vrt">>> ds = gdal.OpenEx(url)>>> image = > >>> ds.ReadAsArray(xoff=0, yoff=0, xsize=64, ysize=64, buf_xsize=10, > >>> buf_ysize=10, resample_alg=gdalconst.GRIORA_NearestNeighbour)>>> > >>> image.shape(3, 10, 10) > > I would be very appreciated if Even could explain why GDAL behaves > differently depending on resample_alg. > > On Thu, Aug 27, 2020 at 6:39 PM Sean Gillies wrote: > >> Hi Denis, Even, >> >> On Thu, Aug 27, 2020 at 8:08 AM Even Rouault >> wrote: >> >>> On jeudi 27 août 2020 15:08:02 CEST Denis Rykov wrote: >>> >>> > I found the culprit. If remove this section from each band definition >>> in >>> >>> > VRT file then everything works fine: >>> >>> > >>> >>> > >>> >>> > dummy.tif >>> >>> > 3 >>> >>> > >> >>> > RasterXSize="40961" RasterYSize="139265" dataType="Byte" /> >>> >>> > >>> >>> > >>> >>> > 0 >>> >>> > 0.0>> >>> >>> >>> I'm perhaps missing something, but the above snippet added by RasterIO >>> doesn't make sense to me. It just sets the pixel at (0,0) to 0. I bet it is >>> completely useless. It would make more sense to have xSize/ySize of DstRect >>> to cover the whole raster. And as is set on the >>> VRTRasterBand, this should be used to fill the target raster anyway. >>> >>> >>> >>> As far as why overviews of the original VRT aren't used is concerned, >>> I'm not sure why. I'd have expected it to work, but I must be missing >>> something. >>> >>> >>> >>> Even >>> >> >> That's a rasterio bug: the background "dummy.tif" source should >> be . I'm working >> on a patch for this. >> >> This VRT background fill is different from nodata. The scale offset could >> be different from the nodata value. BTW, It takes advantage of GDAL >> shortcuts predicated on scale ratio. The "dummy.tif" file is never opened >> and doesn't even exist. >> >> -- >> Sean Gillies >> ___ >> gdal-dev mailing list >> gdal-dev@lists.osgeo.org >> https://lists.osgeo.org/mailman/listinfo/gdal-dev > > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Overviews are not taken into account while reading with specified resampling method
Hi Sean. I patched rasterio as you suggested and intermediate file now looks like this (I'm trying now with a public dataset): https://gist.github.com/drnextgis/3cbdbace7b5b8b80c3c6169b109bf9db But when I read it with GDAL using non-nearest algorithm I'm getting the following error: >>> from osgeo import gdal>>> from osgeo.gdal import gdalconst>>> url = >>> "/tmp/rasterio-boundless.vrt">>> ds = gdal.OpenEx(url)>>> image = >>> ds.ReadAsArray(xoff=0, yoff=0, xsize=64, ysize=64, buf_xsize=10, >>> buf_ysize=10, resample_alg=gdalconst.GRIORA_Cubic) ERROR 4: /tmp/dummy.tif: No such file or directory but with NearestNeighbour it works without error: >>> from osgeo import gdal>>> from osgeo.gdal import gdalconst>>> url = >>> "/tmp/rasterio-boundless.vrt">>> ds = gdal.OpenEx(url)>>> image = >>> ds.ReadAsArray(xoff=0, yoff=0, xsize=64, ysize=64, buf_xsize=10, >>> buf_ysize=10, resample_alg=gdalconst.GRIORA_NearestNeighbour)>>> >>> image.shape(3, 10, 10) I would be very appreciated if Even could explain why GDAL behaves differently depending on resample_alg. On Thu, Aug 27, 2020 at 6:39 PM Sean Gillies wrote: > Hi Denis, Even, > > On Thu, Aug 27, 2020 at 8:08 AM Even Rouault > wrote: > >> On jeudi 27 août 2020 15:08:02 CEST Denis Rykov wrote: >> >> > I found the culprit. If remove this section from each band definition in >> >> > VRT file then everything works fine: >> >> > >> >> > >> >> > dummy.tif >> >> > 3 >> >> > > >> > RasterXSize="40961" RasterYSize="139265" dataType="Byte" /> >> >> > >> >> > >> >> > 0 >> >> > 0.0> >> >> >> I'm perhaps missing something, but the above snippet added by RasterIO >> doesn't make sense to me. It just sets the pixel at (0,0) to 0. I bet it is >> completely useless. It would make more sense to have xSize/ySize of DstRect >> to cover the whole raster. And as is set on the >> VRTRasterBand, this should be used to fill the target raster anyway. >> >> >> >> As far as why overviews of the original VRT aren't used is concerned, I'm >> not sure why. I'd have expected it to work, but I must be missing something. >> >> >> >> Even >> > > That's a rasterio bug: the background "dummy.tif" source should > be . I'm working > on a patch for this. > > This VRT background fill is different from nodata. The scale offset could > be different from the nodata value. BTW, It takes advantage of GDAL > shortcuts predicated on scale ratio. The "dummy.tif" file is never opened > and doesn't even exist. > > -- > Sean Gillies > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Overviews are not taken into account while reading with specified resampling method
On jeudi 27 août 2020 15:08:02 CEST Denis Rykov wrote: > I found the culprit. If remove this section from each band definition in > VRT file then everything works fine: > > > dummy.tif > 3 >RasterXSize="40961" RasterYSize="139265" dataType="Byte" /> > > > 0 > 0.0 is set on the VRTRasterBand, this should be used to fill the target raster anyway. As far as why overviews of the original VRT aren't used is concerned, I'm not sure why. I'd have expected it to work, but I must be missing something. Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Overviews are not taken into account while reading with specified resampling method
I found the culprit. If remove this section from each band definition in
VRT file then everything works fine:
dummy.tif
3
0
0.0https://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 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:
>
>
>
> 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",50],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32647"]]
> 443070.0,1.0,0.0,4366312.0,0.0,-1.0
>
> 0.0
> Red
>
>shared="0">dummy.tif
> 1
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>/>
>/>
> 0
> 0.0
>
>
>shared="0">/vsicurl/https://*.vrt
> 1
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>ySize="139264" />
>yOff="-105176.0" ySize="139264.0" />
> 0.0
>
>
>
> 0.0
> Green
>
>shared="0">dummy.tif
> 2
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>/>
>/>
> 0
> 0.0
>
>
>shared="0">/vsicurl/https://*.vrt
> 2
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>ySize="139264" />
>yOff="-105176.0" ySize="139264.0" />
> 0.0
>
>
>
> 0.0
> Blue
>
>shared="0">dummy.tif
> 3
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>/>
>/>
> 0
> 0.0
>
>
>shared="0">/vsicurl/https://*.vrt
> 3
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>ySize="139264" />
>yOff="-105176.0" ySize="139264.0" />
> 0.0
>
>
>
>
>
>shared="0">/vsicurl/https://*.vrt
> mask,1
>BlockYSize="128" RasterXSize="40961" RasterYSize="139265" dataType="Byte" />
>yOff="0" ySize="139264" />
>yOff="-105176.0" ySize="139264" />
>
>
>
>
> 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
Re: [gdal-dev] Overviews are not taken into account while reading with specified resampling method
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:
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",50],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32647"]]
443070.0,1.0,0.0,4366312.0,0.0,-1.0
0.0
Red
dummy.tif
1
0
0.0
/vsicurl/https://*.vrt
1
0.0
0.0
Green
dummy.tif
2
0
0.0
/vsicurl/https://*.vrt
2
0.0
0.0
Blue
dummy.tif
3
0
0.0
/vsicurl/https://*.vrt
3
0.0
/vsicurl/https://*.vrt
mask,1
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 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
[gdal-dev] Overviews are not taken into account while reading with specified resampling method
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?
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
