Download this notebook from github.


Raster masks from vector masks

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

[1]:
import xarray as xr
import xsar
import numpy as np
import os
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']

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

Adding masks

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

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

Here, we add a ocean_mask dataset variable:

[2]:
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']]

[2]:
<xarray.Dataset>
Dimensions:       (pol: 2, line: 834, sample: 1329)
Coordinates:
  * pol           (pol) object 'VV' 'VH'
  * line          (line) float64 9.5 29.5 49.5 ... 1.663e+04 1.665e+04 1.667e+04
  * sample        (sample) float64 9.5 29.5 49.5 ... 2.655e+04 2.657e+04
    spatial_ref   int64 0
Data variables:
    sigma0_ocean  (pol, line, sample) float64 dask.array<chunksize=(1, 834, 1329), meta=np.ndarray>
    ocean_mask    (line, sample) int8 dask.array<chunksize=(834, 1329), meta=np.ndarray>
Attributes: (12/17)
    name:              SENTINEL1_DS:/tmp/S1B_IW_GRDH_1SDV_20181013T062322_201...
    short_name:        SENTINEL1_DS:S1B_IW_GRDH_1SDV_20181013T062322_20181013...
    product:           GRDH
    safe:              S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_01313...
    swath:             IW
    multidataset:      False
    ...                ...
    footprint:         POLYGON ((-1.371244242212007 48.8181227490355, -4.9521...
    coverage:          168km * 265km (line * sample )
    pixel_line_m:      <xarray.DataArray 'azimuthPixelSpacing' ()>\narray(10....
    pixel_sample_m:    <xarray.DataArray 'groundRangePixelSpacing' ()>\narray...
    orbit_pass:        Descending
    platform_heading:  -165.3526446468305

Masks are available as a shapely object (lon/lat coordinates), with xsar.Sentinel1Meta.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.Sentinel1Meta.ll2coords allow converting lon/lat coordinates on shapely objects.

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

[3]:
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)
[3]: