Windowed Reading with EOReader#
Let’s see how to read a window instead of the 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 performances. Please take it as a hint about relative performances 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", "data", "DS3", "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 than 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/data/DS3/CI/eoreader/others/20201220T104856_L8_200030_OLI_TIRS_window.geojson'
# Here we have the vector directly opened as a geopanas.GeoDataFrame
window_gdf = gpd.read_file(window_path)
window_gdf.to_crs("EPSG:4326").hvplot(coastline="10m")
/opt/conda/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
warnings.warn(f'Downloading: {url}', DownloadWarning)
/opt/conda/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
warnings.warn(f'Downloading: {url}', DownloadWarning)
# 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 construct a Window containing directly the pixels
# It is not necessary to have a window with integers values, they will be rounded when needed.
# However, it is better to control exactly the pixels you want so it is advised to round the values yourself
with rasterio.open(str(prod.get_band_paths([SWIR_1])[SWIR_1])) as ds:
window_pix = from_bounds(*window_bounds, ds.transform)
window_pix
Window(col_off=1932.033839040505, row_off=2154.040731206187, width=3320.4831561867395, height=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 20.9 s, sys: 2.21 s, total: 23.2 s
Wall time: 11.9 s
Time with window#
Use directly the window path.
The keyword used to load with a window is window
.
%time band_arr_win_path = load(prod, window_path)
CPU times: user 2.41 s, sys: 225 ms, total: 2.64 s
Wall time: 1.15 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 2.45 s, sys: 146 ms, total: 2.6 s
Wall time: 1.05 s
%time band_arr_win_bounds = load(prod, window_bounds)
CPU times: user 2.24 s, sys: 213 ms, total: 2.45 s
Wall time: 1.01 s
%time band_arr_win_pix = load(prod, window_pix)
CPU times: user 2.41 s, sys: 132 ms, total: 2.54 s
Wall time: 1.07 s
Some checks#
✅ All windowed bands have the same shape
✅ Windowed bands are smaller that raw band
assert band_arr_win_path.shape == band_arr_win_gdf.shape == band_arr_win_bounds.shape == band_arr_win_pix.shape
band_arr_win_path.shape
(1, 2311, 3321)
band_arr.shape
(1, 7811, 7681)
Plot#
The image with the
[blue - white - yellow]
colomap is the MNDWI band.The image with the
[green - blue]
colomap 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
%%opts QuadMesh [fig_size=500]
# 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])
)