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

Create logger#

# Create logger
import logging
from sertit import logs

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

Open the COSMO-SkyMed 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 COSMO-SkyMed data
path = os.path.join("/home", "data", "DATA", "PRODS", "COSMO", "1st_GEN", "1001512-735097")

# Open your product
prod = Reader().open(path, remove_tmp=True)
prod
eoreader.CskProduct 'CSKS4_DGM_B_HI_09_HH_RA_FF_20201008224018_20201008224025'
Attributes:
	condensed_name: 20201008T224018_CSK_HH_HI_DGM
	path: /home/data/DATA/PRODS/COSMO/1st_GEN/1001512-735097
	constellation: COSMO-SkyMed
	sensor type: SAR
	product type: DGM
	default pixel size: 2.5
	default resolution: 5.0
	acquisition datetime: 2020-10-08T22:40:18.446381
	band mapping:
		HH: HH
		HH_DSPK: HH_DSPK
	needs extraction: True
	orbit direction: ASCENDING
# Plot the quicklook
prod.plot()
/opt/conda/lib/python3.11/site-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
/opt/conda/lib/python3.11/site-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
/opt/conda/lib/python3.11/site-packages/rioxarray/_io.py:1132: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  warnings.warn(str(rio_warning.message), type(rio_warning.message))  # type: ignore
/opt/conda/lib/python3.11/site-packages/rasterio/__init__.py:314: NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific.
  dataset = writer(
../_images/679cc34339883c44ce6e010af99dcc113b0a7c1277cc9e33366ff4dfbfad00d6.png
# 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)
Acquisition datetime: 2020-10-08 22:40:18.446381
Condensed name: 20201008T224018_CSK_HH_HI_DGM
SNAP Release version 9.0.0
SNAP home: /opt/snap/bin/..
SNAP debug: null
SNAP log level: WARNING
Java home: /opt/snap/jre/jre
Java version: 1.8.0_242
Processors: 16
Max memory: 40.9 GB
Cache size: 23.0 GB
Tile parallelism: 14
Tile size: 512 x 512 pixels

To configure your gpt memory usage:
Edit snap/bin/gpt.vmoptions

To configure your gpt cache size and parallelism:
Edit .snap/etc/snap.properties or gpt -c ${cachesize-in-GB}G -q ${parallelism} 
Executing processing graph
Copernicus_DSM_COG_10_N15_00_E108_00_DEM.tif
....10%....20%....30%....40%....50%....60%....70%....80%....90% done.
/opt/conda/lib/python3.11/site-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/snap/s2tbx [reason: /opt/snap/s2tbx]
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/snap/s3tbx [reason: /opt/snap/s3tbx]
2023-11-02 15:51:00,059 - [DEBUG] - Pre-process SAR image
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/snap/s2tbx [reason: /opt/snap/s2tbx]
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/snap/s3tbx [reason: /opt/snap/s3tbx]
2023-11-02 16:07:11,798 - [DEBUG] - Converting DIMAP to GeoTiff
<Axes: >
../_images/58d65a82a3a353cefad140ea6597601e3dab14cac30652460008b12dc3ee124b.png

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", "data", "DS2", "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))
['HH', 'HH_DSPK', 'SLOPE']
# Load those bands as a xarray.Dataset, with a 20m pixel size
band_ds = prod.load(ok_bands, pixel_size=20.)
band_ds[HH]
Executing processing graph
first_line_time metadata value is null
last_line_time metadata value is null
....10%....20%....30%....40%....50%....60%....70%....80%....90% done.
2023-11-02 16:07:30,658 - [DEBUG] - Loading bands ['HH', 'HH_DSPK']
2023-11-02 16:07:30,662 - [DEBUG] - Despeckling HH
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/snap/s2tbx [reason: /opt/snap/s2tbx]
SEVERE: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Failed to scan /opt/snap/s3tbx [reason: /opt/snap/s3tbx]
2023-11-02 16:08:49,827 - [DEBUG] - Loading DEM bands ['SLOPE']
2023-11-02 16:08:49,828 - [DEBUG] - Warping DEM for 20201008T224018_CSK_HH_HI_DGM
2023-11-02 16:08:49,831 - [DEBUG] - Using DEM: /home/data/DS2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt
2023-11-02 16:08:50,653 - [DEBUG] - Computing slope for 20201008T224018_CSK_HH_HI_DGM
<xarray.DataArray <SarBandNames.HH: 'HH'> (band: 1, y: 2473, x: 2688)>
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
Coordinates:
  * x            (x) float64 2.058e+05 2.059e+05 ... 2.596e+05 2.596e+05
  * y            (y) float64 1.746e+06 1.746e+06 ... 1.697e+06 1.697e+06
    spatial_ref  int64 0
  * band         (band) int64 1
Attributes:
    long_name:         HH
    constellation:     COSMO-SkyMed
    constellation_id:  CSK
    product_path:      /home/data/DATA/PRODS/COSMO/1st_GEN/1001512-735097
    product_name:      CSKS4_DGM_B_HI_09_HH_RA_FF_20201008224018_20201008224025
    product_filename:  1001512-735097
    instrument:        SAR-2000
    product_type:      DGM
    acquisition_date:  20201008T224018
    condensed_name:    20201008T224018_CSK_HH_HI_DGM
    orbit_direction:   ASCENDING
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][:, ::10, ::10].plot()
<matplotlib.collections.QuadMesh at 0x7f59f0193450>
../_images/c79889546f9ccbd8ded067364d22560f0a6cba2215abc13565dce7f3602458a6.png

Stack some data#

# You can also stack those bands
stack = prod.stack(ok_bands)
stack
2023-11-02 16:08:53,400 - [DEBUG] - Loading bands ['HH', 'HH_DSPK']
2023-11-02 16:08:53,435 - [DEBUG] - Loading DEM bands ['SLOPE']
2023-11-02 16:08:53,436 - [DEBUG] - Already existing DEM for 20201008T224018_CSK_HH_HI_DGM. Skipping process.
2023-11-02 16:08:53,437 - [DEBUG] - Already existing slope DEM for 20201008T224018_CSK_HH_HI_DGM. Skipping process.
2023-11-02 16:09:14,902 - [DEBUG] - Stacking
/opt/conda/lib/python3.11/site-packages/xarray/core/indexes.py:662: RuntimeWarning: '<' not supported between instances of 'SarBandNames' and 'SarBandNames', sort order is undefined for incomparable objects.
  new_pd_index = pd_indexes[0].append(pd_indexes[1:])
<xarray.DataArray 'HH_HH_DSPK_SLOPE' (bands: 3, y: 19786, x: 21505)>
dask.array<transpose, shape=(3, 19786, 21505), dtype=float32, chunksize=(1, 2048, 2048), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 2.058e+05 2.058e+05 ... 2.596e+05 2.596e+05
  * y            (y) float64 1.746e+06 1.746e+06 ... 1.697e+06 1.697e+06
    spatial_ref  int64 0
  * bands        (bands) object MultiIndex
  * variable     (bands) object SarBandNames.HH ... DemBandNames.SLOPE
  * band         (bands) int64 1 1 1
Attributes:
    long_name:         HH HH_DSPK SLOPE
    constellation:     COSMO-SkyMed
    constellation_id:  CSK
    product_path:      /home/data/DATA/PRODS/COSMO/1st_GEN/1001512-735097
    product_name:      CSKS4_DGM_B_HI_09_HH_RA_FF_20201008224018_20201008224025
    product_filename:  1001512-735097
    instrument:        SAR-2000
    product_type:      DGM
    acquisition_date:  20201008T224018
    condensed_name:    20201008T224018_CSK_HH_HI_DGM
    orbit_direction:   ASCENDING
# 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 in range(nrows):
    stack[i, ::10, ::10].plot(x="x", y="y", ax=axes[i])
../_images/54ce64c83db0bd2aa6fb82ad064abeefe19411769fb655c84ce99028b174022d.png