Dubois surface inversion tutorial¶

This notebook shows how to run the Dubois surface inversion on the UAVSAR multilooked GRD covariance product for winnip_09002_12060_004_120714_L090_CX_02, together with a matching incidence-angle raster.

How the input product was prepared¶

The example data used here were prepared in three steps:

• Download the free UAVSAR multilooked GRD covariance product for winnip_09002_12060_004_120714_L090_CX_02 from the UAVSAR website.
• Import the covariance product into ESA SNAP and export it to NetCDF-BEAM format.
• Download the matching incidence-angle file for the same scene and use the annotation file (.ann) to retrieve the original incidence-angle raster dimensions.

In the notebook below we read the NetCDF-BEAM covariance product with open_netcdf_beam, load the incidence-angle raster, downsample it to match the multilooked grid, and then pass both inputs to dubois_surface_inversion.

Import packages and set directories¶

In [1]:

%load_ext autoreload
%autoreload 2

from pathlib import Path

import numpy as np
import xarray as xr
from dask.diagnostics import ProgressBar

from polsarpro.io import open_netcdf_beam
from polsarpro.physical_inversion import dubois_surface_inversion

# change these paths to match your local files
product_id = "winnip_09002_12060_004_120714_L090_CX_02"

input_uavsar_nc = Path(f"/data/psp/test_files/{product_id}_ML5X5.nc")
incidence_angle_file = Path(f"/data/psp/test_files/{product_id}.inc")
output_dir = Path("/data/psp/res")
output_dir.mkdir(parents=True, exist_ok=True)
file_out = output_dir / 'dubois_surface_inversion_winnip_python.nc'

# Dubois expects GHz; UAVSAR metadata stores the center wavelength in cm.
with xr.open_dataset(input_uavsar_nc) as raw_ds:
    center_wavelength_cm = raw_ds["metadata"].attrs["Annotation:CenterWavelength"]

freq_ghz = 29.9792458 / float(center_wavelength_cm)

# Same parameter values used in the comparison notebook.
thresh1_db = 0.0
thresh2_db = 0.0
calibration_coeff = 1.0

print(f"Derived radar frequency: {freq_ghz:.6f} GHz")

%load_ext autoreload %autoreload 2 from pathlib import Path import numpy as np import xarray as xr from dask.diagnostics import ProgressBar from polsarpro.io import open_netcdf_beam from polsarpro.physical_inversion import dubois_surface_inversion # change these paths to match your local files product_id = "winnip_09002_12060_004_120714_L090_CX_02" input_uavsar_nc = Path(f"/data/psp/test_files/{product_id}_ML5X5.nc") incidence_angle_file = Path(f"/data/psp/test_files/{product_id}.inc") output_dir = Path("/data/psp/res") output_dir.mkdir(parents=True, exist_ok=True) file_out = output_dir / 'dubois_surface_inversion_winnip_python.nc' # Dubois expects GHz; UAVSAR metadata stores the center wavelength in cm. with xr.open_dataset(input_uavsar_nc) as raw_ds: center_wavelength_cm = raw_ds["metadata"].attrs["Annotation:CenterWavelength"] freq_ghz = 29.9792458 / float(center_wavelength_cm) # Same parameter values used in the comparison notebook. thresh1_db = 0.0 thresh2_db = 0.0 calibration_coeff = 1.0 print(f"Derived radar frequency: {freq_ghz:.6f} GHz")

Derived radar frequency: 1.257500 GHz

Load the input product¶

In [ ]:

C3 = open_netcdf_beam(input_uavsar_nc)

C3 = open_netcdf_beam(input_uavsar_nc)

Load and resample the incidence angle¶

We first open the incidence-angle raster at its original dimensions, then downsample it so it matches the 5x5 geocoded image used by the C3 product.

In [3]:

row_dim, col_dim = C3.dims

# use lines from .ann file (infer from auxiliary file)
# hgt.set_rows                         (pixels)        = 3767                   ; HGT Lines
# hgt.set_cols                         (pixels)        = 11636                  ; HGT Samples
nlines = 3767
ncols = 11636

# The incidence-angle file is a little-endian float32 stream in radians.
# Read only the raster-sized prefix and reshape it to the UAVSAR grid.
incidence_angle_values = np.fromfile(
    incidence_angle_file,
    dtype="<f4",
    count=nlines * ncols,
)

# Reduce to the same dimensions as multilooked 5x5 image.
incidence_angle_values[incidence_angle_values == -10000] = np.nan
incidence_angle = xr.DataArray(
    incidence_angle_values.reshape(nlines, ncols),
    dims=(row_dim, col_dim),
    name="incidence_angle",
).coarsen(y=5, x=5, boundary="trim").mean()

row_dim, col_dim = C3.dims # use lines from .ann file (infer from auxiliary file) # hgt.set_rows (pixels) = 3767 ; HGT Lines # hgt.set_cols (pixels) = 11636 ; HGT Samples nlines = 3767 ncols = 11636 # The incidence-angle file is a little-endian float32 stream in radians. # Read only the raster-sized prefix and reshape it to the UAVSAR grid. incidence_angle_values = np.fromfile( incidence_angle_file, dtype="<f4", count=nlines * ncols, ) # Reduce to the same dimensions as multilooked 5x5 image. incidence_angle_values[incidence_angle_values == -10000] = np.nan incidence_angle = xr.DataArray( incidence_angle_values.reshape(nlines, ncols), dims=(row_dim, col_dim), name="incidence_angle", ).coarsen(y=5, x=5, boundary="trim").mean()

Apply the inversion¶

In [4]:

if file_out.exists():
    file_out.unlink()
with ProgressBar():
    dubois_surface_inversion(
        C3,
        incidence_angle=incidence_angle,
        freq_ghz=freq_ghz,
        thresh1=thresh1_db,
        thresh2=thresh2_db,
        calibration_coeff=calibration_coeff,
    ).to_netcdf(file_out)

if file_out.exists(): file_out.unlink() with ProgressBar(): dubois_surface_inversion( C3, incidence_angle=incidence_angle, freq_ghz=freq_ghz, thresh1=thresh1_db, thresh2=thresh2_db, calibration_coeff=calibration_coeff, ).to_netcdf(file_out)

[########################################] | 100% Completed | 1.55 sms

Load the output file¶

In [5]:

res = xr.open_dataset(file_out)

res = xr.open_dataset(file_out)

You can also plot the main Dubois outputs directly:

In [6]:

res.dubois_ks.plot.imshow(cmap='viridis', vmax=0.5, interpolation='bilinear', figsize=(7,7))

res.dubois_ks.plot.imshow(cmap='viridis', vmax=0.5, interpolation='bilinear', figsize=(7,7))

Out[6]:

<matplotlib.image.AxesImage at 0x7f5bc2df7140>

In [7]:

res.dubois_er.plot.imshow(cmap='viridis', vmax=20, interpolation='bilinear', figsize=(7,7))

res.dubois_er.plot.imshow(cmap='viridis', vmax=20, interpolation='bilinear', figsize=(7,7))

Out[7]:

<matplotlib.image.AxesImage at 0x7f5bb41f3c80>

In [8]:

res.dubois_mv.plot.imshow(cmap='viridis', vmax=100, interpolation='bilinear', figsize=(7,7))

res.dubois_mv.plot.imshow(cmap='viridis', vmax=100, interpolation='bilinear', figsize=(7,7))

Out[8]:

<matplotlib.image.AxesImage at 0x7f5ba00fb770>

In [9]:

res.dubois_mask_valid_in_out.plot.imshow(cmap='viridis', interpolation='bilinear', figsize=(7,7))

res.dubois_mask_valid_in_out.plot.imshow(cmap='viridis', interpolation='bilinear', figsize=(7,7))

Out[9]:

<matplotlib.image.AxesImage at 0x7f5ba00d4230>