In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
import bokeh.io
bokeh.io.output_notebook()
import holoviews as hv
hv.extension('bokeh')
import datashader as ds
from holoviews.operation.datashader import datashade,rasterize
from IPython.display import display, HTML
Loading BokehJS ...
In [3]:
import xsar 
from xsarsea import sigma0_detrend
from xsarsea.streaks import streaks_direction
import time

import xarray as xr
import dask.array as da
import numpy as np
from scipy import signal
import sys
import pandas as pd
import numba
from scipy.ndimage import zoom
import time
In [4]:
t1=time.time()
#filename = '/home/alevieux/Documents/Data/S1A_EW_GRDM_1SDV_20160901T202925_20160901T203030_012864_0144F0_C926.SAFE'
#filename = '/home/oarcher/SAFE/S1A_EW_GRDM_1SDV_20160901T202925_20160901T203030_012864_0144F0_C926.SAFE' 
filename = '/home/oarcher/SAFE/S1A_IW_GRDH_1SDV_20170907T103020_20170907T103045_018268_01EB76_992F.SAFE' # Irma

# input dataset, at 100m resolution
sar_ds = xsar.open_dataset(filename,resolution={'atrack': 10, 'xtrack': 10})[0] # full res (10m for IW_GRDH, 40m for EW_GRDM)

# output dataset, at 1km resolution (it's faster to use GDAL rms than to resample sar_ds)
sar_ds_out = xsar.open_dataset(filename,resolution={'atrack': 100, 'xtrack': 100})[0]
# not all variables are needed in the output dataset
sar_ds_out = sar_ds_out.drop(['gamma0', 'gamma0_raw', 'negz' , 'digital_number'])
In [5]:
# detrend sigma0
sar_ds['sigma0_detrend'] = sigma0_detrend(sar_ds.sigma0,sar_ds.incidence)
sar_ds_out['sigma0_detrend'] = sigma0_detrend(sar_ds_out.sigma0,sar_ds_out.incidence)
In [6]:
# compute streaks
streaks_dir = streaks_direction(sar_ds['sigma0_detrend'])

# store streaks in sar_ds_out. (FIXME: interp on degrees need 'nearest' method, need to find a fix for circular)
sar_ds_out['streaks'] = streaks_dir.interp(atrack=sar_ds_out.atrack, xtrack=sar_ds_out.xtrack, method='nearest')

t2 = time.time()
print("elapsed time: %2.1f" % (t2 - t1))

elapsed time: 9.6
In [7]:
rasterize(hv.Image(sar_ds_out['sigma0_detrend'].sel(pol='VV') , rtol=1)).opts(cmap='gray',frame_width=400, frame_height=400) * hv.VectorField( xr.ufuncs.deg2rad(sar_ds_out['streaks'].sel(pol='VV').isel(atrack=slice(None,None,10),xtrack=slice(None,None,10)))).opts(arrow_heads=False)
Out[7]:
In [8]:
# validation with sarwing streaks (estimWindDirSScor_2)

sarwing_ds = xr.open_dataset('/home/oarcher/Téléchargements/s1a-iw-grd-vv-20170907t103020-20170907t103045-018268-01eb76-001_winddir.nc')
hv.VectorField( xr.ufuncs.deg2rad(streaks_dir.sel(pol='VV'))).opts(arrow_heads=False, title='xsar') + hv.VectorField( xr.ufuncs.deg2rad(sarwing_ds['estimWindDirSScor_2'])[::2,::2]).opts(arrow_heads=False,invert_axes=True, title='sarwing')
Out[8]:
In [ ]: