Getting started¶

This tutorial introduces some basic concepts required to use this new version of PolSARpro in python.

It is important to note that the toolbox relies on Xarray. Here we only describe the basics required to run PolSARpro and refer the reader to https://docs.xarray.dev/en/stable/ for more details.

Loading data¶

First we need to import some data. The demo dataset (TODO: host the data somewhere) has been prepared using SNAP and exported in the NetCDF-BEAM. NetCDF is the format recommended by Xarray as they the same data model.

Let's first open the data using the open_netcdf_beam function.

In [ ]:

from polsarpro.io import open_netcdf_beam
from pathlib import Path

# change to your data path
data_path = Path("/data/psp/test_files/SAN_FRANCISCO_ALOS1_slc.nc")
out_path =  Path("/data/psp/res")
S = open_netcdf_beam(data_path)

from polsarpro.io import open_netcdf_beam from pathlib import Path # change to your data path data_path = Path("/data/psp/test_files/SAN_FRANCISCO_ALOS1_slc.nc") out_path = Path("/data/psp/res") S = open_netcdf_beam(data_path)

Note¶

Currently, this function is able to recognize and open the following data types: Sinclair matrix S, covariance matrices C3, C4 and C2 as well as coherency matrices T3 and T4. Preparation of polarimetric data using SNAP is out of the scope of this tutorial and may be found here: https://step.esa.int/docs/tutorials/S1TBX%20Polarimetry%20with%20Radarsat-2%20Tutorial_v2.pdf

Once the data has been prepared in SNAP, it should be exported in NetCDF-BEAM format from the File > Export menu to be opened in PolSARpro.

Inspect data¶

The variable S has been recognized by the reader as a 2x2 Sinclair matrix. It has been imported in a data structure called xarray.Dataset. In a jupyter notebook it may be inspected as:

In [2]:

print(S)

print(S)

<xarray.Dataset> Size: 736MB
Dimensions:  (y: 18432, x: 1248)
Coordinates:
  * y        (y) int64 147kB 0 1 2 3 4 5 ... 18426 18427 18428 18429 18430 18431
  * x        (x) int64 10kB 0 1 2 3 4 5 6 ... 1241 1242 1243 1244 1245 1246 1247
Data variables:
    hh       (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    hv       (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    vv       (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    vh       (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
Attributes:
    poltype:      S
    description:  Scattering matrix

Each polarimetric channel is represented by a data variable and may be accessed using the following dot notation:

In [3]:

# Access to the hh polarization channel
print(S.hh)

# Access to the hh polarization channel print(S.hh)

<xarray.DataArray 'hh' (y: 18432, x: 1248)> Size: 184MB
dask.array<add, shape=(18432, 1248), dtype=complex64, chunksize=(576, 624), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) int64 147kB 0 1 2 3 4 5 ... 18426 18427 18428 18429 18430 18431
  * x        (x) int64 10kB 0 1 2 3 4 5 6 ... 1241 1242 1243 1244 1245 1246 1247

Convert S to T3¶

Each variable is a dask array divided into smaller chunks, enabling parallel processing. Please note that the data is not loaded into memory at this point. This is called lazy loading and allows to execute processing only when the result is required.

More information on lazy operations may be found here:https://docs.xarray.dev/en/latest/internals/internal-design.html#lazy-indexing-classes

To convert the S matrix to other types, we may use utility functions:

In [4]:

from polsarpro.util import S_to_T3

T3 = S_to_T3(S)

from polsarpro.util import S_to_T3 T3 = S_to_T3(S)

Let's now look at this new variable:

In [5]:

print(T3)

print(T3)

<xarray.Dataset> Size: 828MB
Dimensions:  (y: 18432, x: 1248)
Coordinates:
  * y        (y) int64 147kB 0 1 2 3 4 5 ... 18426 18427 18428 18429 18430 18431
  * x        (x) int64 10kB 0 1 2 3 4 5 6 ... 1241 1242 1243 1244 1245 1246 1247
Data variables:
    m11      (y, x) float32 92MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    m22      (y, x) float32 92MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    m33      (y, x) float32 92MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    m12      (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    m13      (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
    m23      (y, x) complex64 184MB dask.array<chunksize=(576, 624), meta=np.ndarray>
Attributes:
    poltype:      T3
    description:  Coherency matrix (3x3)

The output has been automatically converted to a new polarimetric type T3 representing the coherency matrix. Now the elements of the matrix are accessed as:

In [6]:

print(T3.m12)

print(T3.m12)

<xarray.DataArray 'm12' (y: 18432, x: 1248)> Size: 184MB
dask.array<mul, shape=(18432, 1248), dtype=complex64, chunksize=(576, 624), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) int64 147kB 0 1 2 3 4 5 ... 18426 18427 18428 18429 18430 18431
  * x        (x) int64 10kB 0 1 2 3 4 5 6 ... 1241 1242 1243 1244 1245 1246 1247

For storage optimization, PolSARpro takes advantage of the Hermitian structure of the matrix and stores only diagonal and upper elements.

To ensure its algorithms are applied to valid polarimetric types, PolSARpro looks for a poltype "attribute" -- which is a kind of metadata attached to the dataset, and may be accessed as:

In [7]:

T3.poltype

T3.poltype

Out[7]:

'T3'

A human readable description may also be accessed as:

In [8]:

T3.description

T3.description

Out[8]:

'Coherency matrix (3x3)'

Compute the results¶

To compute the actual result and load the data into memory, we may simply use:

In [9]:

T3 = T3.compute()
print(T3)

T3 = T3.compute() print(T3)

<xarray.Dataset> Size: 828MB
Dimensions:  (y: 18432, x: 1248)
Coordinates:
  * y        (y) int64 147kB 0 1 2 3 4 5 ... 18426 18427 18428 18429 18430 18431
  * x        (x) int64 10kB 0 1 2 3 4 5 6 ... 1241 1242 1243 1244 1245 1246 1247
Data variables:
    m11      (y, x) float32 92MB 0.9602 0.03807 0.1286 ... 0.01147 0.003505
    m22      (y, x) float32 92MB 0.01516 0.00855 0.00198 ... 0.001074 0.02919
    m33      (y, x) float32 92MB 0.004146 0.0003961 ... 0.001152 0.0002931
    m12      (y, x) complex64 184MB (-0.054665774+0.10755201j) ... (0.0074339...
    m13      (y, x) complex64 184MB (0.057566743+0.025828078j) ... (0.0008780...
    m23      (y, x) complex64 184MB (-0.00038435453-0.00791823j) ... (0.00087...
Attributes:
    poltype:      T3
    description:  Coherency matrix (3x3)

Crop, convert and save matrices¶

Now, each variable is a numpy array stored in memory.

Lazy loading allows to define complex processing pipelines without storing intermediate variables in memory. It is even possible to write the result to disk without storing the whole data in memory. For instance one may define:

In [ ]:

from polsarpro.util import boxcar
from polsarpro.io import polmat_to_netcdf

# crop the data using the isel accessor, re-chunk to avoid errors when writing zarr
S_crop = S.isel(y=slice(5000, 12000)).chunk("auto")
T3 = S_to_T3(S_crop)

# save to a file with a custom function 
# that preserves chunks and handles complex valued matrices
polmat_to_netcdf(boxcar(T3, 5, 5), out_path / "T3_crop_box_5x5.nc")

from polsarpro.util import boxcar from polsarpro.io import polmat_to_netcdf # crop the data using the isel accessor, re-chunk to avoid errors when writing zarr S_crop = S.isel(y=slice(5000, 12000)).chunk("auto") T3 = S_to_T3(S_crop) # save to a file with a custom function # that preserves chunks and handles complex valued matrices polmat_to_netcdf(boxcar(T3, 5, 5), out_path / "T3_crop_box_5x5.nc")

Apply the H/A/Alpha decomposition¶

Now let's apply the H/A/Alpha decomposition to our S matrix and write the result to a NetCDF file:

In [11]:

from polsarpro.decompositions import h_a_alpha
from dask.diagnostics import ProgressBar

# optional progress bar to monitor the process
with ProgressBar():
    h_a_alpha(S_crop, boxcar_size=[5, 5]).to_netcdf(out_path / "h_a_alpha.nc")

from polsarpro.decompositions import h_a_alpha from dask.diagnostics import ProgressBar # optional progress bar to monitor the process with ProgressBar(): h_a_alpha(S_crop, boxcar_size=[5, 5]).to_netcdf(out_path / "h_a_alpha.nc")

[########################################] | 100% Completed | 32.27 ss

In [12]:

import xarray as xr
# Let's open the file we have written
res = xr.open_dataset(out_path / "h_a_alpha.nc")
print(res)

import xarray as xr # Let's open the file we have written res = xr.open_dataset(out_path / "h_a_alpha.nc") print(res)

<xarray.Dataset> Size: 105MB
Dimensions:     (y: 7000, x: 1248)
Coordinates:
  * y           (y) int64 56kB 5000 5001 5002 5003 ... 11996 11997 11998 11999
  * x           (x) int64 10kB 0 1 2 3 4 5 6 ... 1242 1243 1244 1245 1246 1247
Data variables:
    entropy     (y, x) float32 35MB ...
    anisotropy  (y, x) float32 35MB ...
    alpha       (y, x) float32 35MB ...
Attributes:
    poltype:      h_a_alpha
    description:  Results of the H/A/Alpha decomposition.

As we can see, the new dataset has a specific poltype and stores the output parameters into labeled variables. As previously they can be accessed via the dot notation e.g. res.alpha. We can now plot the outputs either using matplotlib or xarray plotting capabilities:

Visualize the results¶

In [13]:

res.alpha.plot.imshow()

res.alpha.plot.imshow()

Out[13]:

<matplotlib.image.AxesImage at 0x73db8df83d40>

In [14]:

res.entropy.plot.imshow()

res.entropy.plot.imshow()

Out[14]:

<matplotlib.image.AxesImage at 0x73db8c4569c0>

In [15]:

res.anisotropy.plot.imshow()

res.anisotropy.plot.imshow()

Out[15]:

<matplotlib.image.AxesImage at 0x73db8dda6db0>

Additionally, we can plot points in the H/Alpha plane and save the result as an image:

In [16]:

from polsarpro.util import plot_h_alpha_plane
ax, fig = plot_h_alpha_plane(res, bins=400, min_pts=5)
fig.savefig(out_path / "h_a_alpha_plane.png")

from polsarpro.util import plot_h_alpha_plane ax, fig = plot_h_alpha_plane(res, bins=400, min_pts=5) fig.savefig(out_path / "h_a_alpha_plane.png")