Regridding demo datesets using healpix-resample - RIOMAR

Before starting this notebook, Create temporaly small data set using the noteoobk ‘prep_regrid.ipynb’ To install the healpix-resample tool, do

pip install "git+https://github.com/EOPF-DGGS/healpix-resample.git"

import numpy as np
import xarray as xr

# small_ds=xr.open_zarr('/Users/todaka/data/RIOMAR/very_small.zarr')
small_ds = xr.open_zarr("/Users/todaka/data/RIOMAR/small.zarr")

small_ds

small_ds.temp.isel(s_rho=0).hvplot.quadmesh(
    y="nav_lat_rho",
    x="nav_lon_rho",
    ylim=(46.9, 47.5),
)

tmp = (
    small_ds[["temp", "nav_lon_rho", "nav_lat_rho"]]
    .isel(s_rho=0)
    .stack(point=("y_rho", "x_rho"))
)
tmp = tmp.dropna(dim="point").compute()
lon = tmp["nav_lon_rho"].values.astype(np.float64)
lat = tmp["nav_lat_rho"].values.astype(np.float64)
data = tmp.temp.values.astype(np.float64)
data, lon, lat

Here we chose one simple method, but this should be refined with RIOMAR project¶

jmdelouis pablo

from healpix_resample import BilinearResampler

device = "cpu"
level = 13

# nr = BilinearResampler(
#     lon_deg=lon, lat_deg=lat, level=level,
#     device=device, threshold=0.7, ellipsoid="WGS84",
# )
nr = BilinearResampler(
    lon_deg=lon,
    lat_deg=lat,
    level=level,
    device=device,
    threshold=0.5,
    ellipsoid="WGS84",
)
# from healpix_resample import KNeighborsResampler
# nr = KNeighborsResampler(
#     lon_deg=lon, lat_deg=lat, Npt=2,
#     level=level, device=device, ellipsoid="WGS84",
# )
# from healpix_resample import NearestResampler
# level = 12
# nr = NearestResampler(
#     lon_deg=lon, lat_deg=lat, level=level,
#     device=device, ellipsoid="WGS84",
# )

# from healpix_resample import PSFResampler
# level = 13
# nr = PSFResampler(
#     lon_deg=lon, lat_deg=lat, level=level,
#     device=device, threshold=0.7, ellipsoid="WGS84",
# )

# puis:
result = nr.resample(data, lam=0.1)
data_healpix = result.cell_data
# data_healpix_invert = nr.invert(data_healpix)
cell_ids = nr.get_cell_ids()
# hval.shape#,rval.shape
# cell_ids , data_healpix.shape, data_healpix_invert.shape

import xdggs


def make_xdggs_dataarray_from_cell_ids(
    data,
    cell_ids,
    level=15,
    name="da",
):
    """
    Build an xdggs-compatible DataArray:
      - dim/coord: 'cell_ids'
      - values: 1D array (same length as cell_ids)
      - coord attrs required by xdggs
    """
    cell_ids = np.asarray(cell_ids, dtype=np.int64)

    da = xr.DataArray(
        data,
        dims=("cell_ids",),
        coords={"cell_ids": ("cell_ids", cell_ids)},
        name=name,
    )

    da["cell_ids"].attrs.update(
        {
            "grid_name": "healpix",
            "level": int(level),
            "indexing_scheme": "nested",
            "ellipsoid": "WGS84",
        }
    )
    return da.pipe(xdggs.decode).dggs.assign_latlon_coords()


# data_healpix has shape (time, cell_ids)
# — select first timestep for visualization
make_xdggs_dataarray_from_cell_ids(data_healpix[0], cell_ids, level=level).dggs.explore(
    alpha=0.3
)

hp_max = data_healpix.max() - data.max()
hp_min = data_healpix.min() - data.min()
hp_mean = data_healpix.mean() - data.mean()
hp_max, hp_min, hp_mean