import numpy as np
import xsar as sr
import logging
import rasterio
import pickle
import xarray as xr

import matplotlib.pyplot as plt

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

import time

logging.basicConfig()
logger = logging.getLogger('xsar')
logger.setLevel(logging.DEBUG)

Loading BokehJS ...

def display_data(ds_xsar, ds_sarwing, title, polarization):
    diff = (ds_xsar - ds_sarwing).rename("Diff xsar - sarwingsigma0 denoised %s" % polarization)
    std = np.nanstd(diff)
    mean = np.nanmean(diff)
    plots=[]
    r = np.nanpercentile(ds_sarwing,98)
    #r = 0.71
    plots.append(
        rasterize(hv.Image(
            ds_sarwing
            ).opts(cmap='gray',colorbar=True,tools=['hover'],title="sarwing"))
        )
    plots.append(
        rasterize(hv.Image(
            ds_xsar
            ).opts(cmap='gray',colorbar=True,tools=['hover'],title="xsar"))
        )
    plots.append(
        hv.Scatter(
            (ds_xsar,ds_sarwing),
            kdims='xsar', vdims='sarwing'
            ).opts(alpha=0.01,xlim=(0,r),ylim=(0,r),title='std=%.3f bias=%.3f' % (std,mean)) * hv.Path([(0.,0.),(r,r)]).opts(color='black', line_width=2, alpha=0.2)
        )
    return hv.Layout(plots).cols(2).opts(title=title+" "+polarization+" xsar/sarwing")

filename = '/home/alevieux/Documents/Data/S1A_EW_GRDM_1SDV_20160901T202925_20160901T203030_012864_0144F0_C926.SAFE'

polarization = 'VV'

def vent_fictif(Vt):
    # Vt = Vf for Vf =< 10m/s
    Vf1=Vt
    
    # Vt = 0.0009Vf³-0.0173Vf²+1.0831Vf pour 10 < Vf =< 22
    # Find solution of 3rd degre equation using Cardan method
    a=0.0009
    b=-0.0173
    c=1.0831
    d=-Vt
    
    p= c/a-(b*b)/(3.*a*a)
    q=(2.*b*b*b)/(27.*a*a*a)+d/a-(b*c)/(3.*a*a)
    
    # p>0, thus the equation has a unique solution
    D=(q*q)/4.+(p*p*p)/27.

    X=(-q/2.+np.sqrt(D))**(1./3.)-(q/2.+np.sqrt(D))**(1./3.)
    Vf2=X-b/(3.*a)
    
    # Vt = Vf+ Atan(Vf-22)+3.0382 pour Vf >= 22 m/s
    # inversion : fonction affine obtenue par approximation : la pente à
    # l'infini est 1,  
    # l'ordonnée à l'origine bcst est obtenue en moyennant sur 100 valeurs.
    
    # calcul de l'ordonnée à l'origine, en moyennant 
    x=np.arange(100)+50
    bcst=np.mean(np.arctan(x-22.)+3.0382)
    
    Vf3=Vt-bcst

    Vf=(Vt <= 10.)*Vf1 + np.logical_and(Vt > 10. , Vt <= 28.7439)*Vf2 + (Vt > 28.7439)*Vf3
    
    return Vf

def cmodIfr2(wind_speed, wind_dir, inc_angle):
    """
    ; Renvoie une NRCS prédite par modèle IFR2
    ;
    ; INPUTS:
    ;  inc_angle:
    ;  wind_speed:
    ;  wind_dir:
    ;
    ; OPTIONAL INPUTS:
    ;  wind_speed:
    ;
    """

    wind_speed = vent_fictif(wind_speed)
    
    C = np.zeros(26)

    # init CMOD-IFR2 coef
    C[0]  =  0.0
    C[1]  = -2.437597
    C[2]  = -1.5670307
    C[3]  =  0.3708242
    C[4]  = -0.040590
    C[5]  =  0.404678
    C[6]  =  0.188397
    C[7]  = -0.027262
    C[8]  =  0.064650
    C[9]  =  0.054500
    C[10] =  0.086350
    C[11] =  0.055100
    C[12] =  -0.058450
    C[13] = -0.096100
    C[14] = 0.412754
    C[15] =  0.121785
    C[16] =  -0.024333
    C[17] = 0.072163
    C[18] =  -0.062954
    C[19] =  0.015958
    C[20] = -0.069514
    C[21] = -0.062945
    C[22] =  0.035538
    C[23] = 0.023049
    C[24] =  0.074654
    C[25] =  -0.014713
    
    
    T=inc_angle
    wind=wind_speed

    tetai = (T - 36.)/19.
    xSQ = tetai*tetai
    P0 = 1.
    P1 = tetai
    P2 = (3.*xSQ-1.)/2.
    P3 = (5.*xSQ-3.)*tetai/2.
    ALPH = C[1] + C[2] * P1 + C[3] * P2 + C[4]*P3
    BETA = C[5] + C[6] * P1 + C[7] * P2
    ang = wind_dir
    cosi  = np.cos(np.deg2rad(ang))
    cos2i = 2.*cosi*cosi - 1.
    #b1   = self.funcB1(T, wind)
    ############## fonction funcB1 ###################
    tetamin = 18.
    tetamax = 58.
    tetanor = (2.*T - (tetamin+tetamax))/(tetamax-tetamin)
    vmin = 3.
    vmax = 25.
    vitnor = (2.*wind - (vmax+vmin))/(vmax-vmin)
    pv0 = 1.
    pv1 = vitnor
    pv2 = 2*vitnor*pv1 - pv0
    pv3 = 2*vitnor*pv2 - pv1
    pt0 = 1.
    pt1 = tetanor
    pt2 = 2*tetanor*pt1 - pt0
    pt3 = 2*tetanor*pt2 - pt1
    b1 = C[8] + C[9]*pv1 + (C[10]+C[11]*pv1)*pt1 + (C[12]+C[13]*pv1)*pt2
    ##################################################
    #b2   = self.funcB2(T, wind)
    ############## fonction funcB2 ###################
    tetamin = 18.
    tetamax = 58.
    tetanor = (2.*T - (tetamin+tetamax))/(tetamax-tetamin)
    vmin = 3.
    vmax = 25.
    vitnor = (2.*wind - (vmax+vmin))/(vmax-vmin)
    pv0 = 1.
    pv1 = vitnor
    pv2 = 2*vitnor*pv1 - pv0
    pv3 = 2*vitnor*pv2 - pv1
    pt0 = 1.
    pt1 = tetanor
    pt2 = 2*tetanor*pt1 - pt0
    pt3 = 2*tetanor*pt2 - pt1
    result = C[14] + C[15]*pt1 + C[16]*pt2 + \
             (C[17]+C[18]*pt1+C[19]*pt2)*pv1 + (C[20]+C[21]*pt1+C[22]*pt2)*pv2 + \
             (C[23]+C[24]*pt1+C[25]*pt2)*pv3
    b2=result    
    ##################################################
    
    b0 = np.power(10., (ALPH+BETA*np.sqrt(wind)))
    sig  = b0 * ( 1. + b1*cosi + np.tanh(b2)*cos2i )
    return sig

Code

starttime_all = time.time()
image_low_resolution = sr.open_dataset(filename, resolution={'atrack': 100, 'xtrack': 100},pol_dim=True, resampling=rasterio.enums.Resampling.nearest)[0]
image = sr.open_dataset(filename,pol_dim=True)[0]

image['sigma0_denoised'] = np.maximum(image.sigma0 - image.nesz, 0)

image_low_resolution['sigma0_gmf'] = cmodIfr2(10., 45, image_low_resolution.incidence)
image['sigma0_gmf'] = image_low_resolution.sigma0_gmf.interp(atrack=range(0,image.atrack[-1].values,1), xtrack=range(0,image.xtrack[-1].values,1), method='linear')

image.attrs['sigma0_model_mean'] = image_low_resolution.sigma0_gmf.mean()

image['sigma0_detrend'] = np.sqrt(image.sigma0_denoised / (image.sigma0_gmf / image.sigma0_model_mean))

endtime_all = time.time()
logger.debug('temps d exec : %.1fs' % ( endtime_all - starttime_all))

image

DEBUG:xsar.utils:timing load_digital_number : 0.0s. mem: +8.9Mb
DEBUG:xsar.utils:timing load_geometry : 0.5s. mem: +5.4Mb
DEBUG:xsar.utils:timing luts_from_xml_calibration : 0.0s. mem: +3.8Mb
DEBUG:xsar.utils:timing luts_from_xml_calibration : 0.0s. mem: +3.4Mb
DEBUG:xsar.utils:timing luts_from_xml_noise : 0.3s. mem: +31.9Mb
DEBUG:xsar.utils:timing luts_from_xml_noise : 0.3s. mem: +7.1Mb
DEBUG:xsar.utils:timing load_luts : 0.6s. mem: +46.4Mb
DEBUG:xsar.utils:timing luts_from_xml_anotation : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing luts_from_xml_anotation : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing load_luts : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing load_lon_lat : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing _get_subdatasets : 1.4s. mem: +88.1Mb
DEBUG:xsar.utils:timing apply_cf_convention : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing open_dataset : 1.4s. mem: +89.1Mb
DEBUG:xsar.utils:timing load_digital_number : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing load_geometry : 0.5s. mem: +1.1Mb
DEBUG:xsar.utils:timing luts_from_xml_calibration : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing luts_from_xml_calibration : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing luts_from_xml_noise : 0.3s. mem: +29.9Mb
DEBUG:xsar.utils:timing luts_from_xml_noise : 0.3s. mem: +7.1Mb
DEBUG:xsar.utils:timing load_luts : 0.6s. mem: +37.3Mb
DEBUG:xsar.utils:timing luts_from_xml_anotation : 0.0s. mem: +0.3Mb
DEBUG:xsar.utils:timing luts_from_xml_anotation : 0.0s. mem: +0.5Mb
DEBUG:xsar.utils:timing load_luts : 0.0s. mem: +0.8Mb
DEBUG:xsar.utils:timing load_lon_lat : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing _get_subdatasets : 1.3s. mem: +40.2Mb
DEBUG:xsar.utils:timing apply_cf_convention : 0.0s. mem: +0.0Mb
DEBUG:xsar.utils:timing open_dataset : 1.3s. mem: +40.2Mb
DEBUG:xsar:temps d exec : 2.8s

<xarray.Dataset>
Dimensions:          (atrack: 10760, pol: 2, xtrack: 10403)
Coordinates:
  * atrack           (atrack) int64 0 1 2 3 4 ... 10755 10756 10757 10758 10759
  * xtrack           (xtrack) int64 0 1 2 3 4 ... 10398 10399 10400 10401 10402
  * pol              (pol) object 'VV' 'VH'
Data variables:
    digital_number   (pol, atrack, xtrack) uint16 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    time             (atrack) datetime64[ns] dask.array<chunksize=(1000,), meta=np.ndarray>
    longitude        (atrack, xtrack) float32 dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    latitude         (atrack, xtrack) float32 dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    incidence        (atrack, xtrack) float32 dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    elevation        (atrack, xtrack) float32 dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    sigma0_lut       (pol, atrack, xtrack) float32 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    gamma0_lut       (pol, atrack, xtrack) float32 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    noise_lut        (pol, atrack, xtrack) float32 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    sigma0           (pol, atrack, xtrack) float64 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    nesz             (pol, atrack, xtrack) float32 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    gamma0           (pol, atrack, xtrack) float64 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    negz             (pol, atrack, xtrack) float32 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    sigma0_denoised  (pol, atrack, xtrack) float64 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
    sigma0_gmf       (atrack, xtrack) float32 dask.array<chunksize=(10760, 10403), meta=np.ndarray>
    sigma0_detrend   (pol, atrack, xtrack) float64 dask.array<chunksize=(2, 1000, 1000), meta=np.ndarray>
Attributes:
    footprint:          POLYGON ((-40.27692494993099 34.33904879120566, -35.8...
    coverage:           415km * 434km (xtrack * atrack )
    pixel_xtrack_m:     39.9
    pixel_atrack_m:     40.4
    ipf_version:        2.72
    swath_type:         EW
    polarizations:      VV VH
    product_type:       GRD
    mission:            SENTINEL-1
    satellite:          A
    start_date:         2016-09-01 20:29:25.784044
    stop_date:          2016-09-01 20:30:30.065804
    subdataset:         EW
    geometry:                                                                ...
    Conventions:        CF-1.7
    sigma0_model_mean:  <xarray.DataArray 'sigma0_gmf' ()>\ndask.array<mean_a...

Validation

sarwing_sigma0_denoised_full_resolution_vv = pickle.load(open('/home/alevieux/Documents/Debug/validation_ew/sarwing_sigma0_denoised_ew_2016_09_01_vv.pkl', 'rb'), encoding='bytes')
sarwing_sigma0_denoised_full_resolution_vv_dataset = xr.DataArray(sarwing_sigma0_denoised_full_resolution_vv, dims=("atrack", "xtrack"), coords={"atrack":range(0, 6400, 1), "xtrack":range(0,6400,1)})

# slice large tout les 10 point car plantage de jupyter (timeout)
ds_xsar = image.sigma0_denoised.sel(pol='VV', xtrack=slice(0,6399,10), atrack=slice(0,6399,10))
ds_sarwing = sarwing_sigma0_denoised_full_resolution_vv_dataset.sel(xtrack=slice(0,6399,10), atrack=slice(0,6399,10))

display_data(ds_xsar, ds_sarwing, "Validation sigma0 denoised", polarization)

sarwing_sigma0_gmf_full_resolution_vv = pickle.load(open('/home/alevieux/Documents/Debug/validation_ew/sarwing_nrcs_gmf_ew_2016_09_01_vv.pkl', 'rb'), encoding='bytes')
sarwing_sigma0_gmf_full_resolution_vv_dataset = xr.DataArray(sarwing_sigma0_gmf_full_resolution_vv, dims=("atrack", "xtrack"), coords={"atrack":range(0, 6400, 1), "xtrack":range(0,6400,1)})

# slice large tout les 10 point car plantage de jupyter (timeout)
ds_xsar = image.sigma0_gmf.sel(xtrack=slice(0,6399,10), atrack=slice(0,6399,10))
ds_sarwing = sarwing_sigma0_gmf_full_resolution_vv_dataset.sel(xtrack=slice(0,6399,10), atrack=slice(0,6399,10))

display_data(ds_xsar, ds_sarwing, "Validation sigma0 gmf", polarization)

sarwing_sigma0_detrend_full_resolution_vv = pickle.load(open('/home/alevieux/Documents/Debug/validation_ew/sarwing_sigma0_detrend_ew_2016_09_01_vv.pkl', 'rb'), encoding='bytes')
sarwing_sigma0_detrend_full_resolution_vv_dataset = xr.DataArray(sarwing_sigma0_detrend_full_resolution_vv, dims=("atrack", "xtrack"), coords={"atrack":range(0, 6400, 1), "xtrack":range(0,6400,1)})

# slice large tout les 10 point car plantage de jupyter (timeout)
ds_xsar = image.sigma0_detrend.sel(pol='VV', xtrack=slice(0,6399,10), atrack=slice(0,6399,10))
ds_sarwing = sarwing_sigma0_detrend_full_resolution_vv_dataset.sel(xtrack=slice(0,6399,10), atrack=slice(0,6399,10))

display_data(ds_xsar, ds_sarwing, "Validation sigma0 detrend", polarization)

sarwing_sigma0_gmf_all_vv = pickle.load(open('/home/alevieux/Documents/Debug/validation_ew/sarwing_ew_incidence_line1.pkl', 'rb'), encoding='bytes')
sarwing_sigma0_gmf = np.outer(np.ones((10760//1,), dtype=np.float32), sarwing_sigma0_gmf_all_vv)
sarwing_sigma0_gmf_vv_dataset = xr.DataArray(sarwing_sigma0_gmf, dims=("atrack", "xtrack"), coords={"atrack":range(0, 10760, 1), "xtrack":range(0,10403,1)})

# slice large tout les 10 point car plantage de jupyter (timeout)
ds_xsar = image.sigma0_gmf.sel(xtrack=slice(0,image.xtrack.size,100), atrack=slice(0,image.atrack.size,100))
ds_sarwing = sarwing_sigma0_gmf_vv_dataset.sel(xtrack=slice(0,image.xtrack.size,100), atrack=slice(0,image.atrack.size,100))

display_data(ds_xsar, ds_sarwing, "Validation sigma0 gmf gloabl", polarization)

starttime = time.time()
image.sigma0_model_mean.values
endtime = time.time()
logger.debug('sigma0_model_mean %.1fs' % ( endtime - starttime))

starttime = time.time()
image.sigma0.values
endtime = time.time()
logger.debug('sigma0 %.1fs' % ( endtime - starttime))

starttime = time.time()
image.sigma0_gmf.values
endtime = time.time()
logger.debug('sigma0_gmf %.1fs' % ( endtime - starttime))

starttime = time.time()
image.sigma0_detrend.values
endtime = time.time()
logger.debug('sigma0_detrend %.1fs' % ( endtime - starttime))



starttime = time.time()
image.sigma0_model_mean.values
endtime = time.time()
logger.debug('sigma0_model_mean %.1fs' % ( endtime - starttime))

starttime = time.time()
image.sigma0.sel(pol='VV', xtrack=slice(0,6399,1), atrack=slice(0,6399,1)).values
endtime = time.time()
logger.debug('sigma0 %.1fs' % ( endtime - starttime))

starttime = time.time()
image.sigma0_gmf.sel(xtrack=slice(0,6399,1), atrack=slice(0,6399,1)).values
endtime = time.time()
logger.debug('sigma0_gmf %.1fs' % ( endtime - starttime))

starttime = time.time()
image.sigma0_detrend.sel(pol='VV', xtrack=slice(0,6399,1), atrack=slice(0,6399,1)).values
endtime = time.time()
logger.debug('sigma0_detrend %.1fs' % ( endtime - starttime))

DEBUG:xsar:sigma0_model_mean 0.0s
DEBUG:xsar:sigma0 1.9s
DEBUG:xsar:sigma0_gmf 1.6s
DEBUG:xsar:sigma0_detrend 4.9s
DEBUG:xsar:sigma0_model_mean 0.0s
DEBUG:xsar:sigma0 0.7s
DEBUG:xsar:sigma0_gmf 1.6s
DEBUG:xsar:sigma0_detrend 3.0s

sarwing_sigma0_detrend_full_resolution_vv = pickle.load(open('/home/alevieux/Documents/Debug/validation_ew/sarwing_ew_sigma0_detrend_last_tile.pkl', 'rb'), encoding='bytes')
sarwing_sigma0_detrend_full_resolution_vv_dataset = xr.DataArray(sarwing_sigma0_detrend_full_resolution_vv, dims=("atrack", "xtrack"), coords={"atrack":range(6400, 6400+4352, 1), "xtrack":range(6400,6400+3968,1)})

# slice large tout les 10 point car plantage de jupyter (timeout)
ds_xsar = image.sigma0_detrend.sel(pol='VV', xtrack=slice(6400, 10367,10), atrack=slice(6400,10752,10))
ds_sarwing = sarwing_sigma0_detrend_full_resolution_vv_dataset.sel(xtrack=slice(6400,6400+sarwing_sigma0_detrend_full_resolution_vv_dataset.xtrack.size,10), atrack=slice(6400,6400+sarwing_sigma0_detrend_full_resolution_vv_dataset.atrack.size,10))

display_data(ds_xsar, ds_sarwing, "Validation sigma0 detrend last tile", polarization)

rasterize(hv.Image(
            image.sigma0_detrend.sel(pol='VV')
            ).opts(cmap='gray',colorbar=True,tools=['hover'],title="xsar"))

rasterize(hv.Image(
            image.sigma0_denoised.sel(pol='VV')
            ).opts(cmap='gray',colorbar=True,tools=['hover'],title="xsar"))