search

umr-lops / xsar

Synthetic Aperture Radar (SAR) Level-1 GRD python mapper for efficient xarray/dask based processing
https://cyclobs.ifremer.fr/static/sarwing_datarmor/xsar/
MIT License
25 stars 8 forks source link

allow projection from rasterfile to xsar dataset (ecmwf, ice, etc ...) #30

Closed oarcher closed 2 years ago

oarcher commented 2 years ago

We can currently add masks from shapefiles to the dataset with https://github.com/umr-lops/xsar/blob/c204981aa77a80994bfbb034e98d5eb061cf26fb/src/xsar/sentinel1_meta.py#L370-L373

we have also to provide a similar feature for raster files, like netcdf (ecmwf, ice, GEBCO, etc ...).

Each ancillary file should provide a regular lon/lat grid (ie EPSG:4326), or have a CRS from rasterio

If possible, we will also have to define a way to automatically generate the ancillary file name from the SAR date aquisition.

oarcher commented 2 years ago

@vincelhx , can you post here the draft code we made together to map ecmwf file to xsar dataset ?

Put also general overview figure comparison, and local ones that show pixel shift with sarwing

vincelhx commented 2 years ago

Below is a draft code for the add of ecmwf variables (i.e wind speed) into xsar dataset : 

import xsar
import xarray as xr
import rioxarray 
import holoviews as hv
import datetime
hv.extension('bokeh')
filename = '/home/datawork-cersat-public/project/mpc-sentinel1/data/esa/sentinel-1a/L1/IW/S1A_IW_GRDH_1S/2021/252/S1A_IW_GRDH_1SDV_20210909T130650_20210909T130715_039605_04AE83_C34F.SAFE'

sar_ds = xsar.open_dataset(filename, resolution='1000m')
ecwmf_file = '//home/datawork-cersat-public/provider/ecmwf/forecast/hourly/0100deg/netcdf_light/2021/252/ECMWF_FORECAST_0100_202109091300_10U_10V.nc'
ecmwf_ds = xr.open_dataset(ecwmf_file).isel(time=0)
ecmwf_ds.attrs['time'] = datetime.datetime.fromtimestamp(ecmwf_ds.time.item() // 1000000000)
ecmwf_ds = ecmwf_ds.drop('time')
lon = sar_ds.longitude % 360
lat = sar_ds.latitude
U = ecmwf_ds['10U'].interp(Longitude=lon, Latitude=lat)
V = ecmwf_ds['10V'].interp(Longitude=lon, Latitude=lat)
windspeed = np.sqrt(U**2 + V**2)

sar_ds['ecmwf_windspeed'] = windspeed.drop(['Longitude','Latitude'])
sar_ds['ecmwf_windspeed'].attrs.update({k: ecmwf_ds.attrs[k] for k in ['title','institution','time']})
sar_ds['ecmwf_windspeed'].attrs['file'] = ecwmf_file

In order to check if ecmwf_windspeed is in accordance to sarwing_L2 owiEcmwfWindSpeed, we plotted a comparison : 

L2_path  = "/home/datawork-cersat-public/cache/public/ftp/project/sarwing/processings/c39e79a/default/sentinel-1a/L1/IW/S1A_IW_GRDH_1S/2021/252/S1A_IW_GRDH_1SDV_20210909T130650_20210909T130715_039605_04AE83_C34F.SAFE/s1a-iw-owi-xx-20210909t130650-20210909t130715-039605-04AE83.nc"
L2_ds = xr.open_dataset(L2_path)
s1meta = xsar.Sentinel1Meta(filename)

#transform L2 (owiLon,owiLat) coords to (atracks,xtracks)
atracks_2D, xtracks_2D = s1meta.ll2coords(L2_ds['owiLon'],L2_ds['owiLat'])
atracks = np.round(np.mean(atracks_2D, axis=1)).astype(int)
xtracks = np.round(np.mean(xtracks_2D, axis=0)).astype(int)
L2_ds = L2_ds.rename(owiAzSize= "atrack",owiRaSize = "xtrack").assign_coords(atrack = atracks,xtrack = xtracks)
sar_ds["owiEcmwfWindSpeed"] = L2_ds.owiEcmwfWindSpeed.interp(atrack = sar_ds.atrack, xtrack = sar_ds.xtrack,method='nearest')

display : 

hv.Image(sar_ds.owiEcmwfWindSpeed).opts(cmap='jet',colorbar=True, title = "owiEcmwfWindSpeed") +  hv.Image(sar_ds.ecmwf_windspeed).opts(cmap='jet',colorbar=True,title="ecmwf_windspeed") + hv.Image(sar_ds.sigma0.isel(pol=0)).opts(cmap='gray') + hv.Image(sar_ds.owiEcmwfWindSpeed-sar_ds.ecmwf_windspeed).opts(cmap='jet',colorbar=True,title="owiEcmwfWindSpeed - ecmwf_windspeed")

image

display with zoom: 

hv.Image(sar_ds.owiEcmwfWindSpeed).opts(cmap='jet',colorbar=True, title = "owiEcmwfWindSpeed",clim=(0,5)) +  hv.Image(sar_ds.ecmwf_windspeed).opts(cmap='jet',colorbar=True,title="ecmwf_windspeed",clim=(0,5)) + hv.Image(sar_ds.sigma0.isel(pol=0)).opts(cmap='gray') + hv.Image(sar_ds.owiEcmwfWindSpeed-sar_ds.ecmwf_windspeed).opts(cmap='jet',colorbar=True,title="owiEcmwfWindSpeed - ecmwf_windspeed")

image

Comments :

oarcher commented 2 years ago

As a first step, we will close this issue once the following is ok:

However, we will probably need further work:

oarcher commented 2 years ago

image

display with zoom:

image

@vincelhx , can you post here the same validation plot, after you updated your code from #32 ?

vincelhx commented 2 years ago

With @oarcher modifications, the code to load ecmwf variables to xsar dataset is simplified : 

xsar.Sentinel1Meta.set_raster('ecmwf_0100_1h','/home/datawork-cersat-public/provider/ecmwf/forecast/hourly/0100deg/netcdf_light/%Y/%j/ECMWF_FORECAST_0100_%Y%m%d%H%M_10U_10V.nc')

s1meta = xsar.Sentinel1Meta(filename)
sar_ds = xsar.open_dataset(s1meta, resolution='1000m')

leading to 3 new variables in the dataset: ecmwf_0100_1h_U10,ecmwf_0100_1h_V10 and ecmwf_0100_1h_WSPD.

image

Despite the new interpolation method RectBivariateSpline, hatchings are still visible on the difference between sarwing wind speed and wind speed we compute. We think that these hatchings are caused by a difference between our interpolation methods: the current one seems to be a cubic interpolation of an irregular grid to a regular grid.

While zooming on the minimum wind speed we observe that RectBivariateSpline interpolation increased the spatial difference between the 2 supposed eyes

image

On the same area, we used a different colorbar (cmap = "glasbey") to compare the 2 interpolations method around the eye. RectBivariateSpline interpolation makes more sense :

image