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') <https://scitools.org.uk/cartopy/docs/latest/matplotlib/feature_interface.html#cartopy.feature.NaturalEarthFeature>`__

Mask sigma0 over land

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

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

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

<xarray.DataArray 'sigma0_ocean' (pol: 2, atrack: 16685, xtrack: 26584)>
dask.array<transpose, shape=(2, 16685, 26584), dtype=float64, chunksize=(1, 5000, 5000), chunktype=numpy.ndarray>
Coordinates:
  * pol      (pol) object 'VV' 'VH'
  * atrack   (atrack) float64 0.5 1.5 2.5 3.5 ... 1.668e+04 1.668e+04 1.668e+04
  * xtrack   (xtrack) float64 0.5 1.5 2.5 3.5 ... 2.658e+04 2.658e+04 2.658e+04

Adding masks

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

Here, we add a ocean mask:

import cartopy
s1meta = xsar.Sentinel1Meta(xsar.get_test_file('S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_013130_018428_Z000.SAFE'))
s1meta.set_mask_feature('ocean', cartopy.feature.NaturalEarthFeature('physical', 'ocean', '10m'))

Here, we change the default ‘land’ mask for a high resolution shapefile from openstreetmap

s1meta.set_mask_feature('land', os.path.join(xsar.get_test_file('land-polygons-split-4326'),'land_polygons.shp'))

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

s1meta.get_mask('ocean')

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

del ds
ds = xsar.open_dataset(s1meta,resolution={'atrack':10,'xtrack':10})
ds
#ds['ocean_mask']

<xarray.Dataset>
Dimensions:         (atrack: 1668, pol: 2, xtrack: 2658)
Coordinates:
  * pol             (pol) object 'VV' 'VH'
  * atrack          (atrack) float64 5.0 15.0 25.0 ... 1.666e+04 1.668e+04
  * xtrack          (xtrack) float64 5.0 15.0 25.0 ... 2.656e+04 2.658e+04
Data variables: (12/15)
    digital_number  (pol, atrack, xtrack) uint16 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
    time            (atrack) datetime64[ns] 2018-10-13T06:23:22.317102080 ......
    longitude       (atrack, xtrack) float32 dask.array<chunksize=(1668, 2658), meta=np.ndarray>
    latitude        (atrack, xtrack) float32 dask.array<chunksize=(1668, 2658), meta=np.ndarray>
    land_mask       (atrack, xtrack) int8 dask.array<chunksize=(1668, 2658), meta=np.ndarray>
    ocean_mask      (atrack, xtrack) int8 dask.array<chunksize=(1668, 2658), meta=np.ndarray>
    ...              ...
    sigma0_raw      (pol, atrack, xtrack) float64 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
    nesz            (pol, atrack, xtrack) float64 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
    gamma0_raw      (pol, atrack, xtrack) float64 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
    negz            (pol, atrack, xtrack) float64 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
    sigma0          (pol, atrack, xtrack) float64 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
    gamma0          (pol, atrack, xtrack) float64 dask.array<chunksize=(1, 1668, 2658), meta=np.ndarray>
Attributes: (12/15)
    ipf:               2.91
    platform:          SENTINEL-1B
    swath:             IW
    product:           GRDH
    pols:              VV VH
    name:              SENTINEL1_DS:/tmp/S1B_IW_GRDH_1SDV_20181013T062322_201...
    ...                ...
    coverage:          168km * 265km (atrack * xtrack )
    pixel_atrack_m:    100.9310682812434
    pixel_xtrack_m:    99.7062188310889
    orbit_pass:        Descending
    platform_heading:  -165.3526446468305
    Conventions:       CF-1.7

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:

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

land_poly = s1meta.get_mask('land')
land_poly_coords = s1meta.ll2coords(s1meta.get_mask('land'))

ds['sigma0_ocean'] = xr.where(ds['land_mask'], np.nan, ds['sigma0']).transpose(*ds['sigma0'].dims)

rasterize(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)

Performance comparison

# Following code is to compare performances for different methods: it's not usefull for the end user
import pandas as pd
import datetime
index = ['default', 'ne_10m_land', 'gshhs_full_land' , 'osm']
columns = ['feature_f', 'feature_init_time', 'land_footprint_time', 'land_footprint_ax_time', 'rasterize_time']

land_features_df = pd.DataFrame(index=index, columns=columns)
land_features_df.loc['default', 'feature_f'] = lambda: cartopy.feature.LAND
land_features_df.loc['ne_10m_land', 'feature_f'] = lambda: cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')
land_features_df.loc['gshhs_full_land', 'feature_f'] = lambda: cartopy.feature.GSHHSFeature(scale='full')
land_features_df.loc['osm', 'feature_f'] = lambda: os.path.join(xsar.get_test_file('land-polygons-split-4326'),'land_polygons.shp')

raster_shape = s1meta.rio.shape
for feature_idx, feature_rows in land_features_df.iterrows():
    t1 = datetime.datetime.now()
    s1meta.set_mask_feature('land',feature_rows['feature_f']())
    land_features_df.at[feature_idx, 'feature_init_time'] = datetime.datetime.now() - t1
    t1 = datetime.datetime.now()
    land_footprint = s1meta.get_mask('land')
    land_features_df.at[feature_idx, 'land_footprint_time'] = datetime.datetime.now() - t1
    t1 = datetime.datetime.now()
    land_footprint_ax = s1meta.ll2coords(land_footprint)
    land_features_df.at[feature_idx, 'land_footprint_ax_time'] = datetime.datetime.now() - t1
    ds = xsar.open_dataset(s1meta)
    t1 = datetime.datetime.now()
    raster = ds['land_mask'].persist()
    del raster
    land_features_df.at[feature_idx, 'rasterize_time'] = datetime.datetime.now() - t1
land_features_df

	feature_f	feature_init_time	land_footprint_time	land_footprint_ax_time	rasterize_time
default	<function <lambda> at 0x7ff8bd5ee3a0>	0:00:00.014885	0:00:00.033060	0:00:00.004468	0:00:01.538120
ne_10m_land	<function <lambda> at 0x7ff8bd63ab80>	0:00:00.000326	0:00:06.071577	0:00:00.011172	0:00:01.465608
gshhs_full_land	<function <lambda> at 0x7ff8bd63a550>	0:00:00.000546	0:00:47.892499	0:00:00.685888	0:00:03.393890
osm	<function <lambda> at 0x7ff8bd63aa60>	0:00:01.313052	0:00:03.469343	0:00:05.918310	0:00:15.526386