Download this notebook from github.

Metadata vector mask and data raster mask

xsar.SentinelMeta can handle several masks.

A mask is a shapely polygon within the xsar.SentinelMeta.footprint. Default mask is ‘land’, but the user can add customs masks by providing a cartopy.feature.Feature or a shapefile to xsar.SentinelMeta.set_mask_feature.

A mask can be retrieved with xsar.SentinelMeta.get_mask.

When building an xarray dataset from a xsar.SentinelMeta object with xsar.open_dataset, masks are rasterized, to add a variable like ‘land_mask’ to the final dataset. (Not yet implemented)

[1]:

import geopandas as gpd
import pandas as pd
import datetime
import xsar
from osgeo import ogr, gdal
import shapely
import cartopy
import os
import datetime
import rasterio
import rasterio.features

[2]:

# use holoviews for plots
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

Metadata vector mask

xsar.SentinelMeta.footprint is the aquisition footprint, e.g. a rectangle around the acquisition.

Metadata masks are a shapely.geometry object, computed from the intersection between xsar.SentinelMeta.footprint and a cartopy.feature.Feature object.

[3]:

# here is the acquisition footprint, whithout land information
filename = xsar.get_test_file('S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_013130_018428_Z000.SAFE')
s1meta = xsar.SentinelMeta(filename)
s1meta.footprint

[3]:
../_images/examples_mask_4_0.svg

A mask is a shape within xsar.SentinelMeta.footprint that can be retrieved by xsar.SentinelMeta.get_mask

[4]:

default_land_footprint = s1meta.get_mask('land')
default_land_footprint

[4]:
../_images/examples_mask_6_0.svg

A mask can be changed or added with xsar.SentinelMeta.set_mask_feature.

Here, we are defining a new ‘ocean’ mask:

[5]:

s1meta.set_mask_feature('ocean', cartopy.feature.OCEAN)
s1meta.get_mask('ocean')

[5]:
../_images/examples_mask_8_0.svg

High resolution mask

Default cartopy.feature.LAND is low resolution, but xsar.SentinelMeta.set_mask_feature allow to give better resolution features with cartopy.feature.Feature or a user defined shapefile.

  • NaturalEarthFeature

With cartopy.feature.NaturalEarthFeature

[6]:

# the 'land' mask is ovewritten
s1meta.set_mask_feature('land', cartopy.feature.NaturalEarthFeature('physical', 'land', '10m'))
ne_10m_land_footprint = s1meta.get_mask('land')
ne_10m_land_footprint

[6]:
../_images/examples_mask_11_0.svg

  • GSHHSFeature

With cartopy.feature.GSHHSFeature

[7]:

s1meta.set_mask_feature('land',cartopy.feature.GSHHSFeature(scale='full'))
gshhs_land_footprint = s1meta.get_mask('land')
gshhs_land_footprint

[7]:
../_images/examples_mask_13_0.svg

  • Custom feature from shapefile (openstreetmap)

Cartopy has no backends for openstreetmap shapefiles, but we can give any shapefile that can be read by geopandas.read_file to xsar.SentinelMeta.set_mask_feature

Optimisation note

  • Only shapes intersecting with footprint will be read in the shapefile. So performances will be better if shapes are split in smaller shapes (e.g. Eurasian polygon is not an huge polygon)

  • A splitted version from openstreetmap shapefile has been used to speedup access.

[8]:

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


[8]:
../_images/examples_mask_15_0.svg

Resolutions comparison

(Use the mouse to change zoom level)

[9]:

(gv.Shape(default_land_footprint).opts(title='default') \
+ gv.Shape(ne_10m_land_footprint).opts(title='ne_10m') \
+ gv.Shape(gshhs_land_footprint).opts(title='gshhs_full') \
+ gv.Shape(osm_land_footprint).opts(title='osm')).cols(2)

[9]:

Coordinates Reference System

A metadata mask is defined in longitude/latitude reference (EPSG:4326), but the sar data has no CRS. However, sar data has Ground Control Points (GCP), that can be used to convert longitude and latitude to and from atrack and xtrack, with resp s1meta.ll2coords and s1meta.ll2coords

[10]:

osm_land_footprint_ax = s1meta.ll2coords(osm_land_footprint)
osm_land_footprint_ax

[10]:
../_images/examples_mask_19_0.svg

So we are now able to plot the mask over the data

[11]:

ds = xsar.open_dataset(s1meta)
rasterize(hv.Image(((ds.digital_number.sel(pol='VH'))**(1/5))).persist().opts(cmap='gray',colorbar=True,tools=['hover'],title="xsar")) \
* hv.Path(osm_land_footprint_ax).opts(color='lightgreen',width=800,height=800)

[11]:

Data raster mask

Data raster mask is currently not implemented. However, here is a method to get it. ( Warning : this is in-memory, non optimized for dask ! )

[12]:

raster = rasterio.features.rasterize([osm_land_footprint_ax],out_shape=(ds.incidence.shape[1],ds.incidence.shape[0]), all_touched=True).T
ds['land_mask'] = ( {'atrack': raster.shape[1] , 'xtrack': raster.shape[0]}, raster)
ds['dn_masked'] = (ds.digital_number.sel(pol='VH')**(1/5)).where(ds['land_mask'] == 0).persist()
rasterize(hv.Image(ds['dn_masked'])).opts(cmap='gray', width=800,height=800) \
* hv.Path(osm_land_footprint_ax).opts(color='lightgreen',width=800,height=800)

[12]:

Performance comparison

[13]:

# Following code is to compare performances for different methods: it's not usefull for the end user

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
    t1 = datetime.datetime.now()
    raster = rasterio.features.rasterize([land_footprint],out_shape=raster_shape)
    land_features_df.at[feature_idx, 'rasterize_time'] = datetime.datetime.now() - t1
land_features_df


[13]:
feature_f feature_init_time land_footprint_time land_footprint_ax_time rasterize_time
default <function <lambda> at 0x7fef4aaa50d0> 0:00:00.013614 0:00:00.011148 0:00:00.005042 0:00:00.757961
ne_10m_land <function <lambda> at 0x7fef4aaa5280> 0:00:00.000244 0:00:05.648501 0:00:00.009166 0:00:00.715172
gshhs_full_land <function <lambda> at 0x7fef4aaa5310> 0:00:00.000245 0:00:29.468447 0:00:00.653650 0:00:00.814667
osm <function <lambda> at 0x7fef4aaa5670> 0:00:01.223544 0:00:03.556667 0:00:05.825217 0:00:01.019148