SAR example

SAR example#

Let’s use EOReader with SAR data.

Warning: SAR data is processed with SNAP, so be sure to have it installed and that GPT is in your path.

Imports#

import os
import matplotlib.pyplot as plt

# EOReader
from eoreader.reader import Reader
from eoreader.bands import VV, HH, VV_DSPK, HH_DSPK, HILLSHADE, SLOPE, to_str
from eoreader.env_vars import DEM_PATH

/home/docs/checkouts/readthedocs.org/user_builds/eoreader/envs/stable/lib/python3.9/site-packages/dask/dataframe/__init__.py:42: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

Create logger#

# Create logger
import logging
from sertit import logs

logger = logging.getLogger("eoreader")
logs.init_logger(logger)

Open the Capella product#

Please be aware that:

  • EOReader will orthorectify your SAR data to get UTM tiles.

  • complex data is not handled as is, EOReader will convert them to ground range.

# First of all, we need some VHR data, let's use some CAPELLA data
path = os.path.join("/home", "prods", "CAPELLA", "CAPELLA_C02_SS_SLC_HH_20210926061004_20210926061020")

# Open your product
prod = Reader().open(path, remove_tmp=True)
prod

2025-05-05 11:28:13,331 - [WARNING] - There is no existing products in EOReader corresponding to /home/prods/CAPELLA/CAPELLA_C02_SS_SLC_HH_20210926061004_20210926061020.
2025-05-05 11:28:13,333 - [INFO] - Your given path may not be a satellite image. If it is, maybe the product isn't handled by EOReader. If you are sure this product is handled, it is either corrupted or you may need to go deeper in the filetree to find the correct path to give.
2025-05-05 11:28:13,334 - [DEBUG] - Please look at what folder you should give to EOReader by accessing the documentation: https://eoreader.readthedocs.io/latest/main_features.html#recognized-paths

# Plot the quicklook
prod.plot()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[4], line 2
      1 # Plot the quicklook
----> 2 prod.plot()

AttributeError: 'NoneType' object has no attribute 'plot'

# Get the band information
prod.bands

eoreader.SarBand 'HH'
Attributes:
	id: HH
	eoreader_name: HH
	gsd (m): 2.5
	asset_role: intensity
eoreader.SarBand 'HH_DSPK'
Attributes:
	id: HH_DSPK
	eoreader_name: HH_DSPK
	gsd (m): 2.5
	asset_role: intensity

# Print some data
print(f"Acquisition datetime: {prod.datetime}")
print(f"Condensed name: {prod.condensed_name}")

# Open here some more interesting geographical data: extent and footprint
extent = prod.extent()
footprint = prod.footprint()

base = extent.plot(color='cyan', edgecolor='black')
footprint.plot(ax=base, color='blue', edgecolor='black', alpha=0.5)

For SAR data, the footprint needs the orthorectified data ! For that, SNAP uses its own DEM, but you can change it when positionning the EOREADER_SNAP_DEM_NAME environment variable. Available DEMs are:

  • ACE2_5Min

  • ACE30

  • ASTER 1sec GDEM

  • Copernicus 30m Global DEM

  • Copernicus 90m Global DEM

  • GETASSE30

  • SRTM 1Sec HGT

  • SRTM 3Sec

  • External DEM

Warning: If External DEM is set, you must specify the DEM you want by positioning the EOREADER_DEM_PATH to a DEM that can be read by SNAP.

Load bands#

# Set the DEM
os.environ[DEM_PATH] = os.path.join("/home", "ds2_db2", "BASES_DE_DONNEES", "GLOBAL", "COPDEM_30m", "COPDEM_30m.vrt")

# Select some bands you wish to load without knowing if they exist
bands = [VV, HH, VV_DSPK, HH_DSPK, HILLSHADE, SLOPE]

# Only keep those selected
ok_bands = [band for band in bands if prod.has_band(band)]

# This product does not have VV band and HILLSHADE band cannot be computed from SAR band
print(to_str(ok_bands))

# Load those bands as a xarray.Dataset, with a 30 m pixel size
band_ds = prod.load(ok_bands, pixel_size=30.)
band_ds[HH]
Some SAR band (i.e. COSMO) may contain null pixels that are not really nodata (but very low values like water).

This can lead the Terrain Correction step to create large nodata area when projecting on a DEM.

If it happens, you can set the keyword SAR_INTERP_NA to True when loading or stacking SAR data to fill these area with interpolated data.

from eoreader.keywords import SAR_INTERP_NA
band_dict = prod.load(
    ok_bands, 
    pixel_size=20.,
    **{SAR_INTERP_NA: True}
)

# Plot a subsampled version
band_ds[SLOPE][:, ::2, ::2].plot()

Stack some data#

# You can also stack those bands
stack = prod.stack(ok_bands)
stack

# Plot a subsampled version
nrows = len(stack)
fig, axes = plt.subplots(nrows=nrows, figsize=(3 * nrows, 6 * nrows), subplot_kw={"box_aspect": 1})

for i, band_arr in enumerate(stack):
    band_arr[..., ::10, ::10].plot(x="x", y="y", ax=axes[i])
    axes[i].set_title(band_arr.name)
Contents