Download this notebook from github.

Wind speed inversion from level-1 product

Warning

Use of ancillary wind

On this notebook, we changed from ancillary_wind = -np.conj(sarwing_ds.owi_ancillary_wind) to sarwing_ds.owi_ancillary_wind ; then it won’t match sarwing results anymore

[1]:

import xsarsea
from xsarsea import windspeed
import xarray as xr
import numpy as np
import holoviews as hv
hv.extension('bokeh')

import os,sys,re, cv2

[2]:

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

Requirements for inversion

[3]:

import xsar

Getting metadata

[4]:

safe_path = xsarsea.get_test_file("S1A_EW_GRDM_1SDV_20230217T002336_20230217T002412_047268_05AC30_Z005.SAFE")
s1meta = xsar.Sentinel1Meta(safe_path)

[5]:

safe_path

[5]:

'/tmp/S1A_EW_GRDM_1SDV_20230217T002336_20230217T002412_047268_05AC30_Z005.SAFE'

[6]:

from xsarsea.utils import _load_config
config = _load_config()

land mask: not applied in the doc

getting ancillary data

[7]:

s1meta.set_raster('ecmwf_0100_1h',config["path_ecmwf_0100_1h"])
import datetime
for ecmwf_name in ['ecmwf_0100_1h' ]:
    ecmwf_infos = s1meta.rasters.loc[ecmwf_name]
    ecmwf_file = ecmwf_infos['get_function'](ecmwf_infos['resource'], date=datetime.datetime.strptime(s1meta.start_date, '%Y-%m-%d %H:%M:%S.%f'))
    ecmwf_file = xsarsea.get_test_file(ecmwf_file[1].split('/')[-1],iszip=False)
    map_model = { '%s_%s' % (ecmwf_name, uv) : 'model_%s' % uv for uv in ['U10', 'V10'] }

[8]:

map_model

[8]:

{'ecmwf_0100_1h_U10': 'model_U10', 'ecmwf_0100_1h_V10': 'model_V10'}

[9]:

s1meta.rasters.at["ecmwf_0100_1h","resource"] = ecmwf_file

Mapping model & adding ancillary wind

[10]:

### Loading dataset & merging ancillary
xsar_obj = xsar.Sentinel1Dataset(s1meta, resolution='1000m')

[11]:

dataset = xsar_obj.datatree['measurement'].to_dataset()
dataset = dataset.rename(map_model)

## variables of interest

creation of variables of interest for inversion here we could add a land/ice mask.

[12]:

### Variables of interest
#xsar_obj.dataset['land_mask'].values = cv2.dilate(xsar_obj.dataset['land_mask'].values.astype('uint8'),np.ones((3,3),np.uint8),iterations = 3)
#xsar_obj.dataset['sigma0_ocean'] = xr.where(xsar_obj.dataset['land_mask'], np.nan, xsar_obj.dataset['sigma0'].compute()).transpose(*xsar_obj.dataset['sigma0'].dims)
#xsar_obj.dataset['sigma0_ocean'] = xr.where(xsar_obj.dataset['sigma0_ocean'] <= 0, 1e-15, xsar_obj.dataset['sigma0_ocean'])

[13]:

dataset['sigma0_ocean'] = xr.where(dataset['sigma0'] <= 0, 1e-15, xsar_obj.dataset['sigma0'])
dataset['ancillary_wind_direction'] = (90. - np.rad2deg(np.arctan2(dataset.model_V10, dataset.model_U10)) + 180) % 360
dataset['ancillary_wind_speed'] = np.sqrt(dataset['model_U10']**2+dataset['model_V10']**2)
dataset['ancillary_wind'] = dataset.ancillary_wind_speed * np.exp(1j * xsarsea.dir_meteo_to_sample(dataset.ancillary_wind_direction, dataset.ground_heading)) # ref antenna

[14]:

hv.Image(dataset['sigma0_ocean'].sel(pol='VH')).opts(colorbar=True,cmap='binary',width=425, height=400, tools = ['hover'], title = "sigma0 VH")

[14]:

Inversion

inversion parameters

[15]:

apply_flattening = True
GMF_VV_NAME = "gmf_cmod5n"
GMF_VH_NAME = "gmf_s1_v2"

apply flattening or not

[16]:

nesz_cr = dataset.nesz.isel(pol=1) #(no_flattening)
if apply_flattening :
    dataset=dataset.assign(nesz_VH_final=(['line','sample'],windspeed.nesz_flattening(nesz_cr, dataset.incidence)))
    dataset['nesz_VH_final'].attrs["comment"] = 'nesz has been flattened using windspeed.nesz_flattening'
else :
    dataset=dataset.assign(nesz_VH_final=(['line','sample'],nesz_cr.values))
    dataset['nesz_VH_final'].attrs["comment"] = 'nesz has not been flattened'

compute dsig_cr (mix between polarisations) using the last version : “gmf_s1_v2”

[17]:

dsig_cr = windspeed.get_dsig("gmf_s1_v2", dataset.incidence,dataset.sigma0_ocean.sel(pol='VH'),dataset.nesz_VH_final)

retrieve windspeed & direction in dfferent polarizations

define kwargs to parameter LUT resolution

[18]:

kwargs = {"wspd_step" : 0.1, "inc_step" : 0.1, "phi_step" : 1.0, "resolution": "high"}

copol and dual polarization

[19]:

wind_co, wind_dual = windspeed.invert_from_model(
        dataset.incidence,
        dataset.sigma0_ocean.isel(pol=0),
        dataset.sigma0_ocean.isel(pol=1),
        ancillary_wind=dataset['ancillary_wind'],
        dsig_cr = dsig_cr,
        model=(GMF_VV_NAME,GMF_VH_NAME),
        ** kwargs)

dataset["windspeed_co"] = np.abs(wind_co)
dataset["windspeed_dual"] = np.abs(wind_dual)
dataset['winddir_co'] = (90 - np.angle(wind_co, deg=True) + dataset.ground_heading) % 360
dataset['winddir_dual'] = (90 - np.angle(wind_dual, deg=True) + dataset.ground_heading) % 360

cross polarization

[20]:

windspeed_cr = windspeed.invert_from_model(
    dataset.incidence.values,
    dataset.sigma0_ocean.isel(pol=1).values,
    dsig_cr = dsig_cr.values,
    model=GMF_VH_NAME,
    ** kwargs)

dataset = dataset.assign(
    windspeed_cross=(['line', 'sample'], windspeed_cr))

/home/vincelhx/Documents/autoentreprise/IFREMER/libs/fork_xsarsea_vinc/xsarsea/src/xsarsea/windspeed/windspeed.py:80: UserWarning: Unable to check sigma0 pol. Assuming  VH
  warnings.warn(

[21]:

windspeed_co = dataset.windspeed_co.compute()  # compute data if needed
windspeed_cross = dataset.windspeed_cross
windspeed_dual = dataset.windspeed_dual.compute()

windspeeed illustration

[22]:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

projection = ccrs.PlateCarree()
fig, axs = plt.subplots(1, 3, figsize=(18, 6), subplot_kw={'projection': projection})

# Wind Speed Co-Pol
img0 = axs[0].pcolormesh(dataset.longitude, dataset.latitude, windspeed_co, cmap='jet', vmin=0, vmax=80, transform=ccrs.PlateCarree())
axs[0].add_feature(cfeature.COASTLINE)
axs[0].add_feature(cfeature.BORDERS, linestyle=':')
axs[0].set_title('Wind Speed Co-Pol')
gl0 = axs[0].gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')
gl0.top_labels = gl0.right_labels = False

# Wind Speed Cr-Pol
img1 = axs[1].pcolormesh(dataset.longitude, dataset.latitude, windspeed_cross, cmap='jet', vmin=0, vmax=80, transform=ccrs.PlateCarree())
axs[1].add_feature(cfeature.COASTLINE)
axs[1].add_feature(cfeature.BORDERS, linestyle=':')
axs[1].set_title('Wind Speed Cr-Pol')
gl1 = axs[1].gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')
gl1.top_labels = gl1.right_labels = False

# Wind Speed Dual-Pol
img2 = axs[2].pcolormesh(dataset.longitude, dataset.latitude, windspeed_dual, cmap='jet', vmin=0, vmax=80, transform=ccrs.PlateCarree())
axs[2].add_feature(cfeature.COASTLINE)
axs[2].add_feature(cfeature.BORDERS, linestyle=':')
axs[2].set_title('Wind Speed Dual-Pol')
gl2 = axs[2].gridlines(draw_labels=True, color='gray', alpha=0.5, linestyle='--')
gl2.top_labels = gl2.right_labels = False

fig.colorbar(img0, ax=axs, orientation='horizontal', fraction=0.02, pad=0.1, aspect=40)
plt.show()

../_images/examples_windspeed_retrieval_L1_38_0.png

winddir illustration

[23]:

sub_ds = dataset.sel(pol='VV').isel(sample=slice(None, None, 10), line=slice(None, None, 10))

vectorfield = hv.VectorField(
    (
        sub_ds.sample, sub_ds.line,
        xsarsea.dir_meteo_to_sample(sub_ds.winddir_dual,sub_ds.ground_heading),
        sub_ds.windspeed_dual
    )
)

hv.Image(dataset.windspeed_dual, kdims=['sample','line']).opts(title='speed and dir', clim=(0,50), cmap='jet') * vectorfield

[23]:

[24]:

sub_ds = dataset.sel(pol='VV').isel(sample=slice(None, None, 10), line=slice(None, None, 10))

vectorfield = hv.VectorField(
    (
        sub_ds.sample, sub_ds.line,
        xsarsea.dir_meteo_to_sample(sub_ds.ancillary_wind_direction,sub_ds.ground_heading),
        sub_ds.ancillary_wind_speed
    )
)

hv.Image(dataset.ancillary_wind_speed, kdims=['sample','line']).opts(title='ECMWF speed and dir', clim=(0,50), cmap='jet') * vectorfield

[24]:

[ ]: