map_raster ; changes suggestions

[vincelhx]

I'm using hwrf as a raster with hwrf_0015_3h(fname, **kwargs)from raster_readers.py. However map_raster - which is used to map the raster on the SAR grid - is defined to be used only with a L1. It's a problem since anyone would use this function for a L2 product. It's the case for me, i need to map hwrf on the new L2-GRD.

So first suggestion : addapt the code to make map_raster usable for any dataset.

To use it on my L2-GRD dataset I adapted the map_raster function. I pointed an 'issue' using the interpolation of hwrf. It happened on these lines :

for var in raster_ds:
            raster_da = raster_ds[var].chunk(raster_ds[var].shape)
            # upscale in original projection using interpolation
            # in most cases, RectBiVariateSpline give better results, but can't handle Nans
            if np.any(np.isnan(raster_da)):
                upscaled_da = raster_da.interp(x=lons, y=lats)
            else:
                upscaled_da = map_blocks_coords(
                    xr.DataArray(dims=['y', 'x'], coords={'x': lons, 'y': lats}).chunk(1000),
                    RectBivariateSpline(
                        raster_da.y.values,
                        raster_da.x.values,
                        raster_da.values,
                        kx=3, ky=3
                    )
                )
            upscaled_da.name = var
            # interp upscaled_da on sar grid
            mapped_ds_list.append(
                upscaled_da.interp(
                    x=self._dataset.longitude,
                    y=self._dataset.latitude
                ).drop_vars(['x', 'y'])
            )

As written, RectBiVariateSpline can't handle Nans. But the upper condition to not use RectBivariateSpline is not enough for a local model like ecmwf;

Using RectBiVariateSpline : with the NaN problem we got wrong values (on the right) instead of Nans

Using interp :

I suggest to change the code of map_raster to force the use of "interp" for some rasters

here are my inputs :

path_hwrf = "/home/datawork-lops-siam-moawi/PROJECTS/CTROVAGUES/TC_1KM_WRF_MN/outputs/HWRF/SURFACE_CARTESIAN_CNT/2023/03/30/20230330_120000_herman17s_20230330_120000_f000_hwrf_cart_surface_cnt.nc"
kwargs_read = {'date': datetime.datetime(2023, 3, 30, 12, 0)}
raster_ds = hwrf_0015_3h(resource_dec[1],**kwargs_read)
raster_ds['spd'] = np.sqrt(raster_ds.U10**2+raster_ds.V10**2)

path_l2="/home/datawork-cersat-public/cache/public/ftp/project/L2GRD/prod_v4/RCM1_OK2496353_PK2497250_1_SCLNA_20230330_111923_VV_VH_GRD/rcm1--owi-xx-20230330t111910-20230330t112038-_____-_____.nc"
ds = xr.open_dataset(path_l2)
polygon_sar = wkt.loads(ds.attrs['footprint'])

[agrouaze]

@vincelhx I need clarifications, do you suggest modifications in https://github.com/umr-lops/xsar/blob/b0d5951f6125262eda252642c3dbbf085a415c0b/src/xsar/base_dataset.py#L711 ? Your suggestion is to use https://numpy.org/doc/stable/reference/generated/numpy.interp.html instead of https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html ? Which raster products would use interp and which RectBivariateSpline ? What is the point about Level-2 ? xsar is dealing only with Level-1 SAR products (and especially GRD). The NaN in raster products is causing the artifact you described? It seems to me that it is an extrapolation.

[vincelhx]

@agrouaze i was thinking to make map_raster a standalone function since we regularly copy it to use it in another context (LEVEL2)

[agrouaze]

I agree https://github.com/umr-lops/xsar/blob/b0d5951f6125262eda252642c3dbbf085a415c0b/src/xsar/base_dataset.py#L711 could be removed from xsar to be versioned in a separate repository. The new repository would be a dependency for xsar sphinx documentation compilation (especially projections.ipynb)