Windowed Reading with EOReader#

Let’s see how to read a window instead of whole bands of a given product and the potential time-consumption gain.

Warning: The durations shown hereunder may not be representative of your computer’s performance. Please take it as a hint about relative performance between constellations.

Try with Landsat-8#

Let’s open a Landsat-8 OLI collection 2 tile.

# Imports
import os
import geopandas as gpd
import rasterio
from shapely.geometry import box

from rasterio.windows import Window, from_bounds
from eoreader.reader import Reader
from eoreader.bands import NDWI, SWIR_1

import hvplot.pandas
import hvplot.xarray
hvplot.extension('matplotlib')

# Open the product
folder = os.path.join("/home", "ds2_db3", "CI", "eoreader")
path = os.path.join(folder, "optical", "LC08_L1TP_200030_20201220_20210310_02_T1.tar")
reader = Reader()
prod = reader.open(path)
band = NDWI

Windows in EOReader#

There are several ways to read a band with a window in EOReader:

Using a vector stored on disk
Using a geopandas.GeoDataFrame
Using a tuple containing the bounds (in the same CRS as the bands, i.e UTM accessible through prod.crs())
Using directly a window. Note that if you want to pass pixel bounds, you must first create a rasterio.windows.Window

⚠ Even if a vector is passed to the window keyword, only its extent will be used.

# Here we have a window as a path to a vector
window_path = os.path.join(folder, "others", "20201220T104856_L8_200030_OLI_TIRS_window.geojson")
window_path

'/home/ds2_db3/CI/eoreader/others/20201220T104856_L8_200030_OLI_TIRS_window.geojson'

# Here we have the vector directly opened as a geopandas.GeoDataFrame
window_gdf = gpd.read_file(window_path)
window_gdf.to_crs("EPSG:4326").hvplot(coastline="10m")

# Here we have an array of the bounds (tuple or list also works)
window_bounds = window_gdf.bounds.values[0]
window_bounds

array([ 579046.01517122, 4765387.52325064,  678660.50985682,
       4834693.77806381])

# Here we directly construct a Window containing the pixels
# It is not necessary to have a window with integer values, rasterio handles Windows with floating points values.
# However, it is better to control exactly the pixels you want so it is advised to round the values yourself (see below)
tf, *_ = prod.default_transform()
window_pix = from_bounds(*window_bounds, tf)
window_pix

Window(col_off=np.float64(1932.033839040505), row_off=np.float64(2154.040731206187), width=np.float64(3320.4831561867395), height=np.float64(2310.2084937724576))

# Create function to time it easily
def load(prod, window):
    """
    Function that loads the wanted band over a proposed window and cleans the temporary directory.    
    """
    ds = prod.load(
        band,
        window=window
    )
    prod.clean_tmp()
    return ds[band]

Time without window#

Pass None to load without any window

%time band_arr = load(prod, None)

CPU times: user 46.3 s, sys: 3.49 s, total: 49.8 s
Wall time: 28.5 s

/opt/conda/lib/python3.11/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dest = _reproject(

Time with window#

Use the window path directly. The keyword used to load with a window is window.

%time band_arr_win_path = load(prod, window_path)

CPU times: user 8.45 s, sys: 492 ms, total: 8.94 s
Wall time: 6.11 s

Read with a window in different ways#

Here, try to load with:

geopandas.GeoDataFrame
Array of bounds
rasterio.windows.Window

%time band_arr_win_gdf = load(prod, window_gdf)

CPU times: user 8.3 s, sys: 547 ms, total: 8.85 s
Wall time: 5.99 s

%time band_arr_win_bounds = load(prod, window_bounds)

CPU times: user 8.39 s, sys: 516 ms, total: 8.91 s
Wall time: 6.02 s

%time band_arr_win_pix = load(prod, window_pix)

CPU times: user 24.6 s, sys: 1.38 s, total: 26 s
Wall time: 16.4 s

Some checks#

✅ All windowed bands derived from coordinates have the same shape
✅ Windowed bands are smaller than the raw band

assert band_arr_win_path.shape == band_arr_win_gdf.shape == band_arr_win_bounds.shape
assert band_arr_win_path.shape != band_arr_win_pix.shape

print(f"Shape returned from AOI window (coordinates) {band_arr_win_path.shape=}")
print(f"Shape returned from pixel window (pixels in float) {band_arr_win_pix.shape=}")

Shape returned from AOI window (coordinates) band_arr_win_path.shape=(1, 2311, 3321)
Shape returned from pixel window (pixels in float) band_arr_win_pix.shape=(1, 2310, 3320)

band_arr.shape

(1, 7811, 7681)

Plot#

The image with the [blue - white - yellow] colormap is the MNDWI band.
The image with the [green - blue] colormap is the windowed MNDWI band.
The green line corresponds to the window’s bounds
The blue line is the given vector used as a window

# Compute bounds
bounds = band_arr.rio.bounds()
window_bounds_gpd = gpd.GeoDataFrame(
    geometry=[box(*window_bounds)],
    crs=window_gdf.crs
)

# Plot
(band_arr[0, ::10, ::10].hvplot.quadmesh(
    "x", "y", 
    coastline="10m",
    cmap="bwy_r"
) * band_arr_win_path[0, ::10, ::10].hvplot.quadmesh(
    "x", "y", 
    coastline="10m",
    cmap="GnBu"
) * window_bounds_gpd.hvplot(
    facecolor=(0, 0, 0, 0), 
    edgecolor="g", linewidth=4
) * window_gdf.hvplot(
    facecolor=(0, 0, 0, 0), 
    edgecolor="b", linewidth=4, 
    xlim=(bounds[0], bounds[2]), 
    ylim=(bounds[1], bounds[3])
)).opts(fig_size=500)

Windowed Reading with EOReader

Contents