Download this notebook from github.

Streaks analysis

Streaks analysis is done by Koch (20004) algorithm implementation.

[1]:

# import needed modules
import xsar
import xsarsea
import xsarsea.gradients

import xarray as xr
import numpy as np
import scipy
import os
import time

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

[2]:

# optional debug messages
import logging
logging.basicConfig()
logging.getLogger('xsarsea.gradients').setLevel(logging.INFO) # or .setLevel(logging.DEBUG)

[3]:

xsarsea.windspeed.gmfs.GmfModel.activate_gmfs_impl()
xsarsea.windspeed.available_models()

[3]:
alias pol model
gmf_cmod5 cmod5 VV <GmfModel('gmf_cmod5') pol=VV>
gmf_cmod5n cmod5n VV <GmfModel('gmf_cmod5n') pol=VV>
gmf_cmodifr2 cmodifr2 VV <GmfModel('gmf_cmodifr2') pol=VV>
gmf_rs2_v2 rs2_v2 VH <GmfModel('gmf_rs2_v2') pol=VH>
gmf_s1_v2 s1_v2 VH <GmfModel('gmf_s1_v2') pol=VH>
gmf_rcm_noaa rcm_noaa VH <GmfModel('gmf_rcm_noaa') pol=VH> 
[4]:

# open a file a 100m
filename = xsar.get_test_file('S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_018268_01EB76_Z010.SAFE') # irma
#filename = xsar.get_test_file('S1B_IW_GRDH_1SDV_20181013T062322_20181013T062347_013130_018428_Z000.SAFE') # bz
#filename=xsar.get_test_file('S1B_IW_GRDH_1SDV_20211024T051203_20211024T051228_029273_037E47_Z010.SAFE') # nices streaks
#filename=xsar.get_test_file('S1A_IW_GRDH_1SDV_20170720T112706_20170720T112735_017554_01D5C2_Z010.SAFE') # subswath
#filename = '/home/oarcher/SAFE/S1A_EW_GRDM_1SDV_20181009T234410_20181009T234510_024066_02A153_DB9C.SAFE'
sar_ds = xsar.open_dataset(filename,resolution='100m').isel(line=slice(20,None,None),sample=slice(20,None,None)) # isel to skip bad image edge

# add detrended sigma0
sar_ds['sigma0_detrend'] = xsarsea.sigma0_detrend(sar_ds.sigma0, sar_ds.incidence, model = "gmf_cmod5n")
# apply land mask
land_mask = sar_ds['land_mask'].compute()
sar_ds['sigma0_detrend'] = xr.where(land_mask, np.nan, sar_ds['sigma0_detrend']).transpose(*sar_ds['sigma0_detrend'].dims).compute()

DEBUG:xsarsea.windspeed:timing _gmf_function : 0.84s. mem: +15.0Mb

General overview

Gradients direction analysis is done by moving a window over the image. xsarsea.gradients.Gradients allow multiple windows sizes and resolutions.

sar_ds is a IW_GRDH SAFE with a pixel size of 10m at full resolution. So to compute compute gradients with windows size of 16km and 32km, we need to use windows_sizes=[1600,3200]

sar_ds resolution is 100m, so if we want to compute gradients at 100m an 200m, we need to use downscales_factors=[1,2]

[5]:

gradients = xsarsea.gradients.Gradients(sar_ds['sigma0_detrend'], windows_sizes=[1600,3200], downscales_factors=[1,2])
#gradients = xsarsea.gradients.Gradients(sar_ds['sigma0_detrend'], windows_sizes=[400,800], downscales_factors=[1,2])

# get gradients histograms as an xarray dataset
hist = gradients.histogram

# get orthogonals gradients
hist['angles'] = hist['angles'] + np.pi/2

#mean
hist_mean = hist.mean(['downscale_factor','window_size','pol'])

# mean, and smooth
hist_mean_smooth = hist_mean.copy()
hist_mean_smooth['weight'] = xsarsea.gradients.circ_smooth(hist_mean_smooth['weight'])

# smooth only
hist_smooth = hist.copy()
hist_smooth['weight'] = xsarsea.gradients.circ_smooth(hist_smooth['weight'])

# select histogram peak
iangle = hist_mean_smooth['weight'].fillna(0).argmax(dim='angles')
streaks_dir = hist_mean_smooth.angles.isel(angles=iangle)
streaks_weight = hist_mean_smooth['weight'].isel(angles=iangle)
streaks = xr.merge([dict(angle=streaks_dir,weight=streaks_weight)]).drop('angles')


# convert from image convention (rad=0=line) to geographic convention (deg=0=north)
# select needed variables in original dataset, and map them to streaks dataset
streaks_geo = sar_ds[['longitude','latitude','ground_heading']].interp(
    line=streaks.line,
    sample=streaks.sample,
    method='nearest')

streaks_geo['weight'] = streaks['weight']

# convert directions from image convention to geographic convention
streaks_geo['streaks_dir'] =  xsarsea.dir_sample_to_meteo(np.rad2deg(streaks['angle']), streaks_geo['ground_heading'])

streaks_geo = streaks_geo.compute()

# plot. Note that hv.VectorField only accept radians, and 0 is West, so we need to reconvert degrees to radians when calling ...
gv.tile_sources.Wikipedia * gv.VectorField(
    (
        streaks_geo['longitude'],
        streaks_geo['latitude'],
        np.pi/2 -np.deg2rad(streaks_geo['streaks_dir']),
        streaks_geo['weight']
    )
).opts(pivot='mid', arrow_heads=False, tools=['hover'], magnitude='Magnitude')

INFO:xsarsea.gradients:timing histogram : 6.85s. mem: +26.1Mb
/tmp/ipykernel_998005/3506759536.py:43: GeoviewsDeprecationWarning: 'Wikipedia' is deprecated and will be removed in version 1.11, use 'OSM' instead.
  gv.tile_sources.Wikipedia * gv.VectorField(

[5]:

[ ]:

WARNING: hv.VectorField and gv.VectorField don’t use degrees north convention, but radian convention, with 0 = East or right So, to use them with degrees north, you have to convert them to gradients with

np.pi/2 -np.deg2rad(deg_north)

Digging into intermediate computations

streaks_geo

streaks_geo is a xarray.Dataset, with latitude, longitude and streaks_dir (0=deg north) variables.

It has dims ('line', 'sample'), with a spacing corresponding to the first windows size, according to the window step.

[6]:

streaks_geo

[6]:

<xarray.Dataset>
Dimensions:         (line: 11, sample: 16)
Coordinates:
    spatial_ref     int64 0
  * line            (line) float64 204.5 1.804e+03 ... 1.46e+04 1.62e+04
  * sample          (sample) float64 204.5 1.804e+03 ... 2.26e+04 2.42e+04
Data variables:
    longitude       (line, sample) float64 -67.87 -68.02 ... -70.26 -70.41
    latitude        (line, sample) float64 20.71 20.74 20.77 ... 19.65 19.68
    ground_heading  (line, sample) float32 -169.3 -169.3 ... -169.2 -169.2
    weight          (line, sample) float64 0.001487 0.002755 ... 0.001145
    streaks_dir     (line, sample) float64 -80.56 -85.56 ... -203.0 -193.0
Attributes: (12/15)
    name:              SENTINEL1_DS:/tmp/S1A_IW_GRDH_1SDV_20170907T103020_201...
    short_name:        SENTINEL1_DS:S1A_IW_GRDH_1SDV_20170907T103020_20170907...
    product:           GRDH
    safe:              S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_01826...
    swath:             IW
    multidataset:      False
    ...                ...
    start_date:        2017-09-07 10:30:20.936409
    stop_date:         2017-09-07 10:30:45.935264
    footprint:         POLYGON ((-67.84221143971432 20.72564283093837, -70.22...
    coverage:          170km * 251km (line * sample )
    orbit_pass:        Descending
    platform_heading:  -167.7668824808032

streaks

streaks_geo was computed from streaks (also a xarray.Dataset). The main difference is that the angle variable from streaks is in radians, in image convention (ie rad=0 is in sample direction)

[7]:

streaks

[7]:

<xarray.Dataset>
Dimensions:      (line: 11, sample: 16)
Coordinates:
    spatial_ref  int64 0
  * line         (line) float64 204.5 1.804e+03 3.404e+03 ... 1.46e+04 1.62e+04
  * sample       (sample) float64 204.5 1.804e+03 ... 2.26e+04 2.42e+04
Data variables:
    angle        (line, sample) float64 0.02182 0.1091 0.1091 ... 2.16 1.985
    weight       (line, sample) float64 0.001487 0.002755 ... 0.003254 0.001145

Convertion between image convention and geographic convention

see xsarsea.dir_meteo_to_sample and xsarsea.dir_sample_to_meteo

hist_mean

streaks variable was computed from hist_mean_smooth.

The main difference with streaks variable is that we don’t have a single angle, but a histogram of probability for binned angles

[8]:

hist_mean_smooth

[8]:

<xarray.Dataset>
Dimensions:      (line: 11, sample: 16, angles: 72)
Coordinates:
  * line         (line) float64 204.5 1.804e+03 3.404e+03 ... 1.46e+04 1.62e+04
  * sample       (sample) float64 204.5 1.804e+03 ... 2.26e+04 2.42e+04
  * angles       (angles) float64 0.02182 0.06545 0.1091 ... 3.033 3.076 3.12
    spatial_ref  int64 0
Data variables:
    weight       (line, sample, angles) float64 0.001487 0.001487 ... 0.0001425
    used_ratio   (line, sample) float64 0.25 0.5 0.5 0.5 ... 0.2689 0.3732 0.116

Let’s exctract one histogram at an arbitrary position, and plot the histogram.

We can do this with the regular hv.Histogram function, or use xsarsea.gradients.circ_hist, that might be used with hv.Path to plot the histogram as a circular one.

[9]:

hist_at = hist_mean_smooth['weight'].sel(line=5000,sample=12000,method='nearest')
hv.Histogram( (hist_at.angles, hist_at )) + hv.Path(xsarsea.gradients.circ_hist(hist_at))

/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(
/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  layout_plot = gridplot(

[9]:

xsarsea also provide an interactive drawing class xsarsea.gradients.PlotGradients that can be used to draw the circular histogram at mouse tap. (needs a live notebook)

[10]:

# background image for vectorfield
s0 = sar_ds['sigma0_detrend'].sel(pol='VV')
hv_img = rasterize(hv.Image(s0, kdims=['sample', 'line']).opts(cmap='gray',clim=(0,np.nanpercentile(s0,95))))


plot_mean_smooth = xsarsea.gradients.PlotGradients(hist_mean_smooth)

# get vectorfield, with mouse tap activated
hv_vf = plot_mean_smooth.vectorfield(tap=True)

# connect mouse to histogram
hv_hist = plot_mean_smooth.mouse_histogram()

# notebook dynamic output
hv_hist + hv_img * hv_vf

/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(
/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  layout_plot = gridplot(

[10]:

hist_mean_smooth was smoothed. Let’s try hist_smooth

[11]:

plot_smooth = xsarsea.gradients.PlotGradients(hist_smooth)
hv_vf = plot_smooth.vectorfield()
hv_hist = plot_smooth.mouse_histogram()
hv_hist + (hv_img * hv_vf).opts(legend_position='right', frame_width=300)

/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(
/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  layout_plot = gridplot(

[11]:

Using source keyword for mouse_histogram, we can link several histrograms

[12]:

plot_raw = xsarsea.gradients.PlotGradients(hist)
plot_mean = xsarsea.gradients.PlotGradients(hist_mean)
hv_vf = plot_smooth.vectorfield()

hist_smooth_mean = hist_smooth.mean(['downscale_factor','window_size'])
plot_smooth_mean = xsarsea.gradients.PlotGradients(hist_smooth_mean)

gridspace = hv.GridSpace(kdims=['smooth','mean'])
gridspace[(False,False)] = plot_smooth.mouse_histogram(source=plot_raw)
gridspace[(True,False)] = plot_smooth.mouse_histogram()
gridspace[(True,True)] = plot_smooth.mouse_histogram(source=plot_mean_smooth)
gridspace[(False,True)] = plot_smooth.mouse_histogram(source=plot_mean)
#gridspace[(False,True)] = plot_smooth.mouse_histogram(source=plot_smooth_mean)


gridspace.opts(plot_size=(200,200)) + (hv_img * hv_vf).opts(legend_position='right', frame_height=500)

/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:588: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  plot = gridplot(plots[::-1],
/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:647: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  plot = gridplot([r1, r2])
/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_drag' property; using the latest value
  layout_plot = gridplot(
/home/vincelhx/miniconda3/envs/xsar_N3_local/lib/python3.10/site-packages/holoviews/plotting/bokeh/plot.py:987: UserWarning: found multiple competing values for 'toolbar.active_scroll' property; using the latest value
  layout_plot = gridplot(

[12]: