1. Rasterio | Main
  2. Topics

Decimated and windowed read in rasterio

umbertofilippo@...
 

(Posted also on gis.stackexchange.com: https://gis.stackexchange.com/questions/353794/decimated-and-windowed-read-in-rasterio)

In my workflow using rasterio, I'd like to read an overview from a raster and get only a portion of it through a window at the same time.

Is this possible?

I have a pretty complex script so far and I'm trying to check whether the wrong output is due to this not being possible.

Basically, I am asking if it is possible to do something like this:

data = src.read(out_shape=(1,
                           math.ceil(src.height / 64),
                           math.ceil(src.width / 64)),
                window=window)

where data is a decimated read of src to get source raster overview for factor 64, and window has the width and height whose dimensions are relative to the source raster (not the overview).
Sean Gillies
 

Hi,

On Thu, Mar 12, 2020 at 9:30 AM <umbertofilippo@...> wrote:

(Posted also on gis.stackexchange.com: https://gis.stackexchange.com/questions/353794/decimated-and-windowed-read-in-rasterio)

In my workflow using rasterio, I'd like to read an overview from a raster and get only a portion of it through a window at the same time.

Is this possible?

I have a pretty complex script so far and I'm trying to check whether the wrong output is due to this not being possible.

Basically, I am asking if it is possible to do something like this:

data = src.read(out_shape=(1,
                           math.ceil(src.height / 64),
                           math.ceil(src.width / 64)),
                window=window)

where data is a decimated read of src to get source raster overview for factor 64, and window has the width and height whose dimensions are relative to the source raster (not the overview).


the out_shape and window parameters of a Rasterio dataset's read method are very similar to the nBufXSize, nBufYSize, nXOff, nYOff, nXSize, and nYSize parameters of GDAL's RasterIO function. The window parameters refer to a region of the source raster at its full resolution, overviews are disregarded. The out_shape paramter determines the size of the output buffer and GDAL *may* read from overviews when they are available to reduce the cost of resampling.

Rasterio's read(), like GDAL's RasterIO, abstracts the details of full resolution vs overviews. Depending on the math, you may or may not get exactly the overview level you are looking for. If you absolutely need a particular overview level, I suggest passing an overview_level keyword argument to rasterio.open().

--
Sean Gillies
umbertofilippo@...
 

On Thu, Mar 12, 2020 at 05:29 PM, Sean Gillies wrote:
Depending on the math, you may or may not get exactly the overview level you are looking for. If you absolutely need a particular overview level, I suggest passing an overview_level keyword argument to rasterio.open().
This is very interesting, I could not find the reference to overview_level in the documentation.
Definitely going to test it!
My aim is to read a specific overview (the greatest level) from a raster and resampling it by implementing a SUMMING algorithm.
Although this might sound like I am not supposed to use windows, this process is part of a more complex program where the user can also request the original resolution, and thus for very big dataset reading with windows is much better.
Will come back with my findings, meanwhile can you please point to a documentation if present where the open method is used alng with the overview_level paramenter?

Thank you so much,

Umberto Minora
umbertofilippo@...
 

The overview_level did the trick!

Thank you so much, I was after this for a long time.

Umberto Minora
umbertofilippo@...
 
Edited

Thanks to @Sean Gillies.

Posting here the solution I also posted on gis.stackexchange.com (https://gis.stackexchange.com/a/353869/9518).

Maybe it would be useful for others in search of a solution for a similar problem.

Basically, it is not possible to simultaneously make a windowed AND decimated read, so the answer to my question is NO.

However, there is an undocumented (AFAIK) parameter called overview_level which can be passed to rasterio.open().

This parameter takes the 0-based index of an overview level. So, assuming overviews are present in a source raster, one can do

src = rasterio.open(raster, overview_level=0)

to create a Rasterio dataset from the first level of overview of source raster.

Then, simply do a windowed read in the common format:

for ji, src_window in src.block_windows(1):
    arr = src.read(1, window=src_window)

1 - 5 of 5