Download this notebook from github.

[1]:

import xarray as xr
import xsar
import numpy as np
import os

Raster masks from vector masks

Sentinel 1 example

xsar.Sentinel1Dataset.dataset has a land_mask variable by default, rasterized from cartopy.feature.NaturalEarthFeature(‘physical’, ‘land’, ‘10m’)

Mask sigma0 over land

Replacing sigma0 values with nan over ocean can be done with xarray.where

[2]:

ds = xsar.open_dataset(xsar.get_test_file('S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_013130_018428_Z010.SAFE'), resolution='1000m')

# apply land_mask. The final transpose is done to preserve dimensions ordering
ds['sigma0_land'] = xr.where(ds['land_mask'], ds['sigma0'], np.nan).transpose(*ds['sigma0'].dims)
ds['sigma0_land']

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/numpy/_core/numeric.py:452: RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')

[2]:

<xarray.DataArray 'sigma0_land' (pol: 2, line: 166, sample: 265)> Size: 704kB
dask.array<transpose, shape=(2, 166, 265), dtype=float64, chunksize=(1, 166, 265), chunktype=numpy.ndarray>
Coordinates:
  * line         (line) float64 1kB 49.5 149.5 249.5 ... 1.645e+04 1.655e+04
  * sample       (sample) float64 2kB 49.5 149.5 249.5 ... 2.635e+04 2.645e+04
  * pol          (pol) object 16B 'VV' 'VH'
    spatial_ref  int64 8B 0

Adding masks

Masks can be added with xsar.BaseMeta.set_mask_feature, providing a shapefile or a cartopy.feature.Feature object.

For a default mask for all SAFE, use classmethod xsar.BaseMeta.set_mask_feature

Here, we add a ocean_mask dataset variable:

[3]:

xsar.Sentinel1Meta.set_mask_feature('ocean', os.path.join(xsar.get_test_file('water-polygons-split-4326'), 'water_polygons.shp'))


s1meta = xsar.Sentinel1Meta(xsar.get_test_file('S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_013130_018428_Z010.SAFE'))
ds = xsar.open_dataset(s1meta, resolution='200m')
ds['sigma0_ocean'] = xr.where(ds['ocean_mask'], ds['sigma0'], np.nan).transpose(*ds['sigma0'].dims)
ds[['sigma0_ocean','ocean_mask']]

/home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/testgeoviewsnumpy/lib/python3.11/site-packages/numpy/_core/numeric.py:452: RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')

[3]:

<xarray.Dataset> Size: 19MB
Dimensions:       (pol: 2, line: 834, sample: 1329)
Coordinates:
  * line          (line) float64 7kB 9.5 29.5 49.5 ... 1.665e+04 1.667e+04
  * sample        (sample) float64 11kB 9.5 29.5 49.5 ... 2.655e+04 2.657e+04
  * pol           (pol) object 16B 'VV' 'VH'
    spatial_ref   int64 8B 0
Data variables:
    sigma0_ocean  (pol, line, sample) float64 18MB dask.array<chunksize=(1, 834, 1329), meta=np.ndarray>
    ocean_mask    (line, sample) int8 1MB dask.array<chunksize=(834, 1329), meta=np.ndarray>
Attributes: (12/15)
    name:              SENTINEL1_DS:/home1/scratch/agrouaze/xsardatasync/xsar...
    short_name:        SENTINEL1_DS:S1B_IW_GRDH_1SDV_20181013T062322_20181013...
    product:           GRDH
    safe:              S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_01313...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2018-10-13 06:23:22.317102
    stop_date:         2018-10-13 06:23:47.315227
    footprint:         POLYGON ((-1.371244242212007 48.8181227490355, -4.9521...
    coverage:          168km * 265km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -165.3526446468305

Masks are available as a shapely object (lon/lat coordinates), with xsar.BaseMeta.get_mask:

When using this s1meta object with xsar.Sentinel1Meta.footprint, mask are rasterized over the dataset, and variables postfixed with _mask are created.

Convert a mask to dataset coordinates

xsar.BaseMeta.ll2coords allow converting lon/lat coordinates on shapely objects.

So it’s possible to plot the vector mask on a raster variable:

[4]:

import holoviews as hv
import geoviews.feature as gf
hv.extension('bokeh')
from holoviews.operation.datashader import datashade,rasterize


# get the shapely polygon of the mask (on the footprint, lon/lat)
land_poly = s1meta.get_mask('ocean')

# convert lon/lat to line/sample
land_poly_coords = s1meta.ll2coords(s1meta.get_mask('ocean'))

# plot masked sigma0, with coastline

(
    hv.Image(ds['sigma0_ocean'].sel(pol='VV')).opts(cmap='gray',clim=(0,0.2), colorbar=True,tools=['hover'],title="xsar") \
    * hv.Path(land_poly_coords).opts(color='lightgreen',width=800,height=800)
).opts(invert_axes=True)

/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)

[4]:

RadarSat 2 example

xsar.RadarSat2Dataset.dataset has a land_mask variable by default, rasterized from cartopy.feature.NaturalEarthFeature(‘physical’, ‘land’, ‘10m’)

Mask sigma0 over land

Replacing sigma0 values with nan over ocean can be done with xarray.where

[5]:

ds = xsar.open_dataset(xsar.get_test_file('RS2_OK135107_PK1187782_DK1151894_SCWA_20220407_182127_VV_VH_SGF'), resolution='1000m')

# apply land_mask. The final transpose is done to preserve dimensions ordering
ds['sigma0_land'] = xr.where(ds['land_mask'], ds['sigma0'], np.nan).transpose(*ds['sigma0'].dims)
ds['sigma0_land']

[5]:

<xarray.DataArray 'sigma0_land' (pol: 2, line: 513, sample: 530)> Size: 4MB
dask.array<transpose, shape=(2, 513, 530), dtype=float64, chunksize=(1, 513, 530), chunktype=numpy.ndarray>
Coordinates:
  * line     (line) float64 4kB 9.5 29.5 49.5 ... 1.021e+04 1.023e+04 1.025e+04
  * pol      (pol) <U2 16B 'VV' 'VH'
  * sample   (sample) float64 4kB 9.5 29.5 49.5 ... 1.057e+04 1.059e+04

[6]:

ds['sigma0_land'].sel(pol='VH').plot(vmin=0)

[6]:

<matplotlib.collections.QuadMesh at 0x7b22a2c26110>

RCM example

xsar.RcmDataset.dataset has a land_mask variable by default, rasterized from cartopy.feature.NaturalEarthFeature(‘physical’, ‘land’, ‘10m’)

[7]:

import warnings
warnings.filterwarnings('ignore')

Mask sigma0 over land

Replacing sigma0 values with nan over ocean can be done with xarray.where

[8]:

ds = xsar.open_dataset(xsar.get_test_file('RCM1_OK1050603_PK1050605_1_SC50MB_20200214_115905_HH_HV_Z010'), resolution='1000m')

# apply land_mask. The final transpose is done to preserve dimensions ordering
ds['sigma0_land'] = xr.where(ds['land_mask'], ds['sigma0'], np.nan).transpose(*ds['sigma0'].dims)
ds['sigma0_land']

[8]:

<xarray.DataArray 'sigma0_land' (pol: 2, line: 382, sample: 360)> Size: 2MB
dask.array<transpose, 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
  * pol      (pol) <U2 16B 'HH' 'HV'
  * sample   (sample) float64 3kB 24.5 74.5 124.5 ... 1.792e+04 1.797e+04

[9]:

ds['sigma0_land'].sel(pol='HV').plot(vmin=0)

[9]:

<matplotlib.collections.QuadMesh at 0x7b22b0fbfe50>