Assume longitude is cyclic in `_interp_nominal_lon`

[JoranAngevaare]

trafficstars

A potential solution to #295. Instead of assuming there are a cyclic 360 datapoints (in longitude), assume that the set is cyclic with however many datapoins in there are in the longitude coordinate.

Example

Setup

# Basic imports
import numpy as np
import xarray as xr
from xmip.preprocessing import promote_empty_dims, replace_x_y_nominal_lat_lon, rename_cmip6, correct_coordinates, broadcast_lonlat
import matplotlib.pyplot as plt
import sys
import numpy, xarray, xmip
message = f'{sys.platform}\n{sys.version}'
for b in 'numpy xarray xmip'.split():
    message += f'\n\t{b}=={locals().get(b).__version__}'
print(message)

# Totally arbitrary data
lon=np.linspace(0,360, 513)[:-1]
lat=np.linspace(-90,90, 181)[:-1]
time=np.arange(1)
data=np.arange(len(lat)*len(lon)*len(time)).reshape(len(time), len(lat), len(lon))*lon

# Add some NaN values just as an example
data[:, :, len(lon)//2+30: len(lon)//2+50] = np.nan

ds_dummy = xr.Dataset(
    data_vars = dict(var=(('time','lat', 'lon'), data,)),
    coords = dict(time=time, lon=lon, lat=lat,),
    attrs=dict(source_id='bla'))


def setup_plt():
    """use cartopy to make a "map" for the dummy data. If cartopy is not installed, do nothing"""
    try:
        import cartopy.crs as ccrs
    except (ImportError, ModuleNotFoundError):
        return
    plt.gcf().add_subplot(projection=ccrs.PlateCarree(central_longitude=0.0,))
    ax = plt.gca()

    ax.coastlines()
    ax.gridlines(draw_labels=True)

def recast(data_set):
    ds = rename_cmip6(data_set)
    ds = promote_empty_dims(ds)
    ds = broadcast_lonlat(ds)
    ds = replace_x_y_nominal_lat_lon(ds)
    return ds

which gives

linux
3.8.16 (default, Mar  2 2023, 03:21:46) 
[GCC 11.2.0]
	numpy==1.24.3
	xarray==2023.1.0
	xmip==0.7.1

Starting point:

setup_plt()
ds_dummy['var'].isel(time=0).plot(cbar_kwargs=dict(orientation='horizontal'))

Broken example:

setup_plt()
ds_dum_recast = recast(ds_dummy)
replace_x_y_nominal_lat_lon(ds_dum_recast)['var'].isel(time=0).plot(cbar_kwargs=dict(orientation='horizontal'),)

Patch:

def _interp_nominal_lon_new(lon_1d):
    print('Using altered version')
    x = np.arange(len(lon_1d))
    idx = np.isnan(lon_1d)
    # TODO assume that longitudes are cyclic (i.e., don't)
    ret = np.interp(x, x[~idx], lon_1d[~idx], period=len(lon_1d))
    return ret

xmip.preprocessing._interp_nominal_lon = _interp_nominal_lon_new
setup_plt()
ds_dum_recast = recast(ds_dummy)
ds_dum_recast['var'].isel(time=0).plot(cbar_kwargs=dict(orientation='horizontal', extend='both'),)

[jbusecke]

Reopening this. Again sorry for the long wait (being lone maintainer on a repo is never ideal 😆). This looks like a nice fix @JoranAngevaare , could you add a few test cases (you could take them directly from #295!) in this PR? Then we can def get this merged.

[jbusecke]

@JoranAngevaare again, big thanks and sorry for the delay. Please add yourself in #358. I really think you made significant contributions!