Download this notebook from github.

Advanced explanation for RCM

Contrary to Sentinel-1, Radarsat Constellation Mission doesn’t have the notion of multi dataset

[1]:

import xsar
import geoviews as gv
import holoviews as hv
import geoviews.feature as gf
hv.extension('bokeh')
path = xsar.get_test_file('RCM1_OK1050603_PK1050605_1_SC50MB_20200214_115905_HH_HV_Z010')

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/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)

[2]:

import warnings
warnings.filterwarnings('ignore')

Open a dataset

[9]:

# Define the resolution to load the dataset at a lower resolution (if not specified or None, the dataset is loaded at high resolution)
resolution = '1000m'

# Instanciate a RadarSatDataset object
rcmds = xsar.RcmDataset(dataset_id=path, resolution=resolution)

Get the Dataset object

[10]:

rcmds

[10]:

<RcmDataset full coverage object>

Access the metadata object from the Dataset object

[11]:

rcmds.sar_meta

[11]:

<BlockingActorProxy: RcmMeta>

Access the dataset

In this dataset, we can find variables like latitude, longitude, look up tables (before and after denoising), incidence…

[12]:

rcmds.dataset

[12]:

<xarray.Dataset> Size: 23MB
Dimensions:          (line: 382, pol: 2, sample: 360)
Coordinates:
  * line             (line) float64 3kB 24.5 74.5 124.5 ... 1.902e+04 1.907e+04
  * pol              (pol) <U2 16B 'HH' 'HV'
  * sample           (sample) float64 3kB 24.5 74.5 ... 1.792e+04 1.797e+04
Data variables: (12/21)
    digital_number   (pol, line, sample) uint16 550kB dask.array<chunksize=(1, 382, 360), meta=np.ndarray>
    lines_flipped    bool 1B False
    samples_flipped  bool 1B True
    sampleSpacing    float64 8B 1e+03
    lineSpacing      float64 8B 1e+03
    sigma0_raw       (pol, line, sample) float64 2MB dask.array<chunksize=(1, 382, 360), meta=np.ndarray>
    ...               ...
    ground_heading   (line, sample) float32 550kB dask.array<chunksize=(382, 360), meta=np.ndarray>
    latitude         (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
    longitude        (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
    altitude         (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
    incidence        (line, sample) float32 550kB dask.array<chunksize=(382, 360), meta=np.ndarray>
    elevation        (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
Attributes: (12/18)
    name:                    RCM_DS:/home1/scratch/agrouaze/xsardatasync/xsar...
    short_name:              RCM_DS:RCM1_OK1050603_PK1050605_1_SC50MB_2020021...
    product:                 Z010
    safe:                    RCM1_OK1050603_PK1050605_1_SC50MB_20200214_11590...
    swath:                   SC50MB
    multidataset:            False
    ...                      ...
    stop_date:               2020-02-14 12:00:02.000000
    footprint:               POLYGON ((-84.63397773377159 48.71534850877886, ...
    coverage:                381km * 359km (line * sample )
    pixel_line_m:            20.0
    pixel_sample_m:          20.0
    approx_transform:        |-0.00,-0.00,-84.71|\n|-0.00, 0.00, 48.73|\n| 0....

Variables lines_flippedand samples_flipped are added to the dataset to know if these have been flipped (in order to follow xsar convention).

See RCM Image Product Format Definition (4.2.1) for more explication about the different cases when lines or samples are flipped

Alternatives solutions to open dataset and datatree

[13]:

# Open dataset
xsar.open_dataset(path, resolution=resolution)

[13]:

<xarray.Dataset> Size: 23MB
Dimensions:          (line: 382, pol: 2, sample: 360)
Coordinates:
  * line             (line) float64 3kB 24.5 74.5 124.5 ... 1.902e+04 1.907e+04
  * pol              (pol) <U2 16B 'HH' 'HV'
  * sample           (sample) float64 3kB 24.5 74.5 ... 1.792e+04 1.797e+04
Data variables: (12/21)
    digital_number   (pol, line, sample) uint16 550kB dask.array<chunksize=(1, 382, 360), meta=np.ndarray>
    lines_flipped    bool 1B False
    samples_flipped  bool 1B True
    sampleSpacing    float64 8B 1e+03
    lineSpacing      float64 8B 1e+03
    sigma0_raw       (pol, line, sample) float64 2MB dask.array<chunksize=(1, 382, 360), meta=np.ndarray>
    ...               ...
    ground_heading   (line, sample) float32 550kB dask.array<chunksize=(382, 360), meta=np.ndarray>
    latitude         (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
    longitude        (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
    altitude         (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
    incidence        (line, sample) float32 550kB dask.array<chunksize=(382, 360), meta=np.ndarray>
    elevation        (line, sample) float64 1MB dask.array<chunksize=(382, 360), meta=np.ndarray>
Attributes: (12/18)
    name:                    RCM_DS:/home1/scratch/agrouaze/xsardatasync/xsar...
    short_name:              RCM_DS:RCM1_OK1050603_PK1050605_1_SC50MB_2020021...
    product:                 Z010
    safe:                    RCM1_OK1050603_PK1050605_1_SC50MB_20200214_11590...
    swath:                   SC50MB
    multidataset:            False
    ...                      ...
    stop_date:               2020-02-14 12:00:02.000000
    footprint:               POLYGON ((-84.63397773377159 48.71534850877886, ...
    coverage:                381km * 359km (line * sample )
    pixel_line_m:            20.0
    pixel_sample_m:          20.0
    approx_transform:        |-0.00,-0.00,-84.71|\n|-0.00, 0.00, 48.73|\n| 0....

[14]:

# Open datatree
xsar.open_datatree(path, resolution=resolution)

[14]:

<xarray.DatasetView> Size: 0B
Dimensions:  ()
Data variables:
    *empty*
Attributes: (12/18)
    name:                    RCM_DS:/home1/scratch/agrouaze/xsardatasync/xsar...
    short_name:              RCM_DS:RCM1_OK1050603_PK1050605_1_SC50MB_2020021...
    product:                 Z010
    safe:                    RCM1_OK1050603_PK1050605_1_SC50MB_20200214_11590...
    swath:                   SC50MB
    multidataset:            False
    ...                      ...
    stop_date:               2020-02-14 12:00:02.000000
    footprint:               POLYGON ((-84.63397773377159 48.71534850877886, ...
    coverage:                381km * 359km (line * sample )
    pixel_line_m:            20.0
    pixel_sample_m:          20.0
    approx_transform:        |-0.00,-0.00,-84.71|\n|-0.00, 0.00, 48.73|\n| 0....

How to apply calibration?

All the operation below are already performed by default for GRD products. So, what is following is a simple explanation about how is made calibration.

Load digital numbers

load_digital_number is a function that allows to load digital numbers from tiff files at chosen resolution and return it as a DataArray. Resampling is made thanks to rasterio.open(tiff_file).read. For dual pol products, there is a tiff file for each pol. So that digital numbers are computed for both pol. Posting of lines and samples is computed thanks to Affine.translation and Affine.scale.

[15]:

import rasterio

#Define resampling method (here it is the root mean square from rasterio)
resampling = rasterio.enums.Resampling.rms

#Define the chunks size for line and samples
chunks = {'line': 5000, 'sample': 5000}

[16]:

dn_low_res = rcmds.load_digital_number(resolution=resolution, resampling=resampling, chunks=chunks)
dn_low_res

[16]:

<xarray.Dataset> Size: 556kB
Dimensions:         (line: 382, sample: 360, pol: 2)
Coordinates:
  * line            (line) float64 3kB 24.5 74.5 124.5 ... 1.902e+04 1.907e+04
  * sample          (sample) float64 3kB 24.5 74.5 124.5 ... 1.792e+04 1.797e+04
  * pol             (pol) <U2 16B 'HH' 'HV'
Data variables:
    digital_number  (pol, line, sample) uint16 550kB dask.array<chunksize=(1, 382, 360), meta=np.ndarray>

Get the raw normalized radar cross section

First the calibration look up tables are loaded at the good resolution in xsar.RcmDataset._luts thanks to the method lazy_load_luts. The resampling is made thanks to an interpolation with the method ̀xsar.RcmDataset._interpolate_var.

This last method uses RectBivariateSplinefor the interpolation, but it is necessary that data is expressed as a 2D vector. Here, calibration/noise look up tables and incidence are expressed as 1D vector (range sample). Consequently, we need to convert these in 2D (adding an azimuth dimension dependency) before applying the interpolation. Conversion is made thanks to numpy.tile, using the low resolution lines expressed in the geolocation grid part of the reader; reducing the calculation. A template of a DataArray that uses the posting of digital numbers (with applied resolution) is given on this interpolation function so the result is now at the right resolution.

_apply_calibration_lut is a method that applies a calibration look up table to digital numbers to compute gamma0, beta0 and sigma0 (depending on the variable name passed in argument) and return the result as a DataArray. It applies the following formula :

\[\frac{(digitalNumbers^2)+offset}{Gain}\]

Reference : RCM Image Product Format Definition (7.5.1)

Different resampling method were tried such as scipy.ndimage.gaussian_filter1d that had the convenience to accept 1d vectors. Data was computed with this function and the chosen posting was this of digital numbers. But in order to be homogenous with other operations made in xsar, we chose to keep the solution with RectBivariateSpline.

[17]:

sigma0_raw = rcmds._apply_calibration_lut('sigma0').sigma0_raw
sigma0_raw

[17]:

<xarray.DataArray 'sigma0_raw' (pol: 2, line: 382, sample: 360)> Size: 2MB
dask.array<where, shape=(2, 382, 360), dtype=float64, chunksize=(1, 382, 360), chunktype=numpy.ndarray>
Coordinates:
  * line     (line) float64 3kB 24.5 74.5 124.5 ... 1.902e+04 1.907e+04
  * sample   (sample) float64 3kB 24.5 74.5 124.5 ... 1.792e+04 1.797e+04
  * pol      (pol) <U2 16B 'HH' 'HV'
Attributes:
    pixelFirstLutValue:  18041
    stepSize:            -23
    numberOfValues:      786
    offset:              0.0

[18]:

import matplotlib.pyplot as plt

[19]:

plt.figure(figsize=(13, 6))

plt.subplot(1, 2, 1)
sigma0_raw_cross = sigma0_raw.sel(pol='HV')
plt.pcolor(sigma0_raw_cross, vmax=0.02, cmap='Greys_r')
plt.title('sigma0_raw HV')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

plt.subplot(1, 2, 2)
sigma0_raw_co = sigma0_raw.sel(pol='HH')
plt.pcolor(sigma0_raw_co, vmax=0.7, cmap='Greys_r')
plt.title('sigma0_raw HH')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

[19]:

<matplotlib.colorbar.Colorbar at 0x74c23fde6f10>

How to apply denoising ?

All the operation below are already performed by default for GRD products. So, what is following is a simple explanation about how is made denoising.

How to get the Noise Equivalent Sigma Zero ?

The noise look up tables are loaded at the good resolution in xsar.RcmDataset._noise_luts thanks to the method lazy_load_noise_luts. The resampling is made thanks to an interpolation with the method ̀xsar.RcmDataset._interpolate_var.

Notes: - noise luts are already calibrated so we don’t have to apply a calibration on these - interpolate_var method explained above converts the NoiseLevelValues in linear beacause these are expressed in dB in the reader. The formula used is :

\[NoiseLevelValues_{linear} = 10^\frac{NoiseLevelValues_{dB}}{10}\]

[20]:

nesz_low_res = rcmds.lazy_load_noise_luts().sigma0
nesz_low_res

[20]:

<xarray.DataArray 'sigma0' (pol: 2, line: 382, sample: 360)> Size: 1MB
dask.array<concatenate, shape=(2, 382, 360), dtype=float32, chunksize=(1, 382, 360), chunktype=numpy.ndarray>
Coordinates:
  * line                  (line) float64 3kB 24.5 74.5 ... 1.902e+04 1.907e+04
  * sample                (sample) float64 3kB 24.5 74.5 ... 1.792e+04 1.797e+04
    pixelFirstNoiseValue  int64 8B 18041
    stepSize              int64 8B -23
  * pol                   (pol) <U2 16B 'HH' 'HV'

[21]:

plt.pcolor(nesz_low_res.sel(pol='HH'), cmap='Greys_r')
plt.title('nesz (pol = HH)')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

[21]:

<matplotlib.colorbar.Colorbar at 0x74c23fa18490>

How to get the noise substracted Sigma0

Right now we only have to substract the noise_lut to the raw normalized radar cross section. It is made with the function _add_denoised, that add the variables to the RcmDataset.dataset

[22]:

sigma0 = sigma0_raw - nesz_low_res
sigma0

[22]:

<xarray.DataArray (pol: 2, line: 382, sample: 360)> Size: 2MB
dask.array<sub, shape=(2, 382, 360), dtype=float64, chunksize=(1, 382, 360), chunktype=numpy.ndarray>
Coordinates:
  * line                  (line) float64 3kB 24.5 74.5 ... 1.902e+04 1.907e+04
  * sample                (sample) float64 3kB 24.5 74.5 ... 1.792e+04 1.797e+04
  * pol                   (pol) <U2 16B 'HH' 'HV'
    pixelFirstNoiseValue  int64 8B 18041
    stepSize              int64 8B -23

Comparison between noised sigma0 and noised substracted sigma0

[23]:

plt.figure(figsize=(26, 12))

sigma0_cross = sigma0.sel(pol='HV')
sigma0_co = sigma0.sel(pol='HH')

plt.subplot(2,2,1)
plt.pcolor(sigma0_raw_cross, vmax=0.02, cmap='Greys_r')
plt.title('Sigma0 HV with noise')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

plt.subplot(2,2,3)
plt.pcolor(sigma0_cross, vmax=0.02, cmap='Greys_r')
plt.title('Sigma0 HV without noise')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

plt.subplot(2,2,2)
plt.pcolor(sigma0_raw_co, vmax=0.7, cmap='Greys_r')
plt.title('Sigma0 HH with noise')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

plt.subplot(2,2,4)
plt.pcolor(sigma0_co, vmax=0.7, cmap='Greys_r')
plt.title('Sigma0 HH without noise')
plt.xlabel('samples')
plt.ylabel('lines')
plt.colorbar()

[23]:

<matplotlib.colorbar.Colorbar at 0x74c21fbdbf50>

How to get the incidence ?

Such as noise/calibration look up tables, incidence depends on samples, which are expressed in the reader thanks to step, firstPixelValue and nbOfValues. Samples are computed in the meta class directly on the datatree thanks to the method assign_index.

_load_incidence_from_lut is a function that applies an interpolation with _interpolate_var and then return the incidence at the good resolution

[24]:

incidence = rcmds._load_incidence_from_lut()
incidence

[24]:

<xarray.DataArray "Delayed('RectBivariateSpline-9a031c3c-ab4f-421f-b5-a2d8ad4798bd5df4af75643ee52ef871" (
                                                                                                         line: 382,
                                                                                                         sample: 360)> Size: 550kB
dask.array<Delayed('RectBivariateSpline-9a031c3c-ab4f-421f-b5, shape=(382, 360), dtype=float32, chunksize=(382, 360), chunktype=numpy.ndarray>
Coordinates:
  * line     (line) float64 3kB 24.5 74.5 124.5 ... 1.902e+04 1.907e+04
  * sample   (sample) float64 3kB 24.5 74.5 124.5 ... 1.792e+04 1.797e+04

How to get the elevation ?

To get the incidence, we apply a formula :

\[\theta = \arcsin{[\sin{(Incidence)} . \frac{r}{r + h}]}\]

\[( r \text{ is the earth radius} , h \text{ is the orbit altitude} )\]

2 variables give orbit altitude so we considered theSatelliteHeight (and not thealtitude).

RcmDataset._load_elevation_from_lut permit to calculate the elevation (in degrees).

[25]:

elevation = rcmds._load_elevation_from_lut()
elevation

[25]:

<xarray.DataArray (line: 382, sample: 360)> Size: 1MB
dask.array<degrees, shape=(382, 360), dtype=float64, chunksize=(382, 360), chunktype=numpy.ndarray>
Coordinates:
  * line     (line) float64 3kB 24.5 74.5 124.5 ... 1.902e+04 1.907e+04
  * sample   (sample) float64 3kB 24.5 74.5 124.5 ... 1.792e+04 1.797e+04