Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Assume longitude is cyclic in _interp_nominal_lon #296

Merged
merged 6 commits into from
Jun 25, 2024

Conversation

JoranAngevaare
Copy link
Contributor

@JoranAngevaare JoranAngevaare commented May 30, 2023

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'))

image

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'),)

image

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'),)

image

@JoranAngevaare JoranAngevaare changed the title Assume longitude is cyclic Assume longitude is cyclic in _interp_nominal_lon May 30, 2023
@jbusecke jbusecke reopened this Jan 25, 2024
@jbusecke
Copy link
Owner

jbusecke commented Jan 25, 2024

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
Copy link
Owner

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

@jbusecke jbusecke merged commit a6248b9 into jbusecke:main Jun 25, 2024
7 of 8 checks passed
JoranAngevaare added a commit to JoranAngevaare/xMIP that referenced this pull request Jun 26, 2024
@JoranAngevaare JoranAngevaare mentioned this pull request Jun 26, 2024
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants