[COG] Reading from the overviews in cog

[Even Rouault]

On mardi 12 juin 2018 16:47:08 CEST Grégory Bataille wrote:
> Hey Even,
> 
> Yes I guess ;)
> 
> However, what you say is exactly what I did and what I had hoped but I
> found that it was reading only the first 512x512 pixels of my big raster
> and dumping at full res rather than downsampling it as I expected!

You probably misused the (xsize, ysize) parameters instead of the (buf_xsize, buf_ysize).

(xoff, yoff, xsize, ysize) are in full resolution image coordinates
buf_xisze, buf_ysize is the dimension of the output buffer

> I expected so much the behavior you mention so much that it took me time to
> understand what was actually happening!
> 
> Maybe I did something wrong then. I’ll retry again

Little demo. Let's generate a 1024x1024 raster, with 512x512 overview

$ gdal_translate byte.tif test.tif -outsize 1024 1024
$ gdaladdo test.tif 2

1) Reading whole raster at full resolution:

$ CPL_DEBUG=ON python -c "from osgeo import gdal; ds = gdal.Open('test.tif'); ds.ReadAsArray()"
[...]
GDAL: GDALOpen(test.tif, this=0x1685430) succeeds as GTiff.
[...]

2) Reading whole raster at overview resolution:

$ CPL_DEBUG=ON python -c "from osgeo import gdal; ds = gdal.Open('test.tif'); ds.ReadAsArray(buf_xsize=512, buf_ysize=512)"
[...]
GDAL: GDALOpen(test.tif, this=0xdae480) succeeds as GTiff.
[...]
GTiff: ScanDirectories()
GTiff: Opened 512x512 overview.
[...]

> 
> Thanks for the quick answer as always
> Greg
> 
> On Tue, 12 Jun 2018 at 16:42, Even Rouault <even.rouault at spatialys.com>
> 
> wrote:
> > On mardi 12 juin 2018 16:20:47 CEST Grégory Bataille wrote:
> > > Hello everyone,
> > > 
> > > I’m trying to migrate my code to COG more and more rather than
> > > pre-computing some views/data.
> > > 
> > > I have had some great success with extracting partial full res raster
> > > values and it’s great.
> > > 
> > > One thing I have not been able to do so far though is « reading from the
> > > overview »
> > > Say for example I want to render a 512x512 tiles of data but really
> > 
> > zoomed
> > 
> > > out. Let’s say for example that this tile contains my entire
> > 
> > 10’000x10’000
> > 
> > > pixels geotiff. It’s obviously super inefficient to read from the raw
> > > raster rather than reading from one of the pre-computed overviews.
> > > Is there an API (python) to ask for raster data (between 2 points) (like
> > > ReadAsArray which is the one I currently use) but passing it the output
> > > size so that gdal can figure out which overview it’s better to read
> > > from?
> > 
> > Grégory,
> > 
> > This is a GDAL question in fact :-)
> > 
> > Yes, the full signature of Dataset.ReadAsArray() is
> > 
> >     def ReadAsArray(self, xoff=0, yoff=0, xsize=None, ysize=None,
> >     
> >                     buf_obj=None,
> >                     buf_xsize=None, buf_ysize=None, buf_type=None,
> >                     resample_alg=gdalconst.GRIORA_NearestNeighbour,
> >                     callback=None,
> >                     callback_data=None,
> > 
> >                     interleave='band'):
> > So specify buf_xsize and buf_ysize (or pass a buf_obj of the appropriate
> > output size)
> > 
> > Even
> > 
> > --
> > Spatialys - Geospatial professional services
> > http://www.spatialys.com
> > 
> > --
> 
> Greg


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com