Handling slicing with circular longitude coordinates in xarray

Visualization
1

I’m working on dividing into basins (Pacific, Atlantic, Indian Southern Ocean) and taking some zonal averages in an xarray dataset. My longitude coordinates are currently 0:360E and so slicing the Atlantic is difficult because lon=slice(298, 22) does not make sense without xarray knowing that this is a circular axis. I know that I could switch the coordinates to -180:180E, but then slicing the Pacific becomes an issue.

I came across this discussion of the xr.roll method which sounds like it might be a way forward, but am having trouble implementing it. I would prefer not to end up with duplicate datasets with different coordinates because I will be adding more derived variables throughout the analysis. Does anyone have an example they’d be willing to share? Thanks!

3 Likes
2

Noting that, as far as I understand it from @TomNicholas, the ongoing Explicit Indexes Refactor of xarray aims to make circular indexing easier when complete. Unfortunately I don’t know how to accomplish it today. Tom, maybe you have a suggestion?

2 Likes
3

What errors are you getting?

This is how I roll coordinates in my dataset (from 0-360 to -180-180) and haven’t had any issues yet. Without having to duplicate the dataset, perhaps this can be done once for the Atlantic and slice Pacific and Indian ocean separately?

    ds = ds.assign_coords(longitude=(((ds.longitude + 180) % 360) - 180))
    ds = ds.roll(longitude=int(len(ds['longitude']) / 2), roll_coords=True)
1 Like
4

Thanks for the example code @josephyang ! That works really well, although it still means I have to have two datasets (or must roll and unroll each time I make a new plot). I’ll stay tuned for the update mentioned above.

1 Like
5

I’m curious does pangeo (or anybody other than the IRI DL) recognize and do the “right thing” with the periodic longitude attribute? It’s pretty handy and seems like the right thing. (In annual cycle calculations the Jan-Dec month grid is periodic too)

grid: /X (degree_east) periodic (0) to (2W) by 2.0 N= 180 pts :grid
http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.ERSST/.version5/.anom/

grid: /T (months since 01-Jan) periodic (Jan) to (Dec) by 1.0 N= 12 pts :grid
http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.ERSST/.version5/.anom/yearly-climatology/

1 Like
6

@josephyang This worked well for me. But how would you go back the other direction? I now have some different datasets that are -180 to 180 but I need to get them to 0-360 so that I can slice across the longitudes 146-292E. Thanks!

7

FYI, this functionality will very soon be supported in xarray.

cc @benbovy

2 Likes
8

This should do it:

ds = ds.assign_coords(longitude=((360 + (ds.longitude % 360)) % 360))
ds = ds.roll(longitude=int(len(ds['longitude']) / 2),roll_coords=True)
3 Likes
9

Yes this looks like another great use case that the flexible indexes refactoring should make easier, e.g., via the implementation of some kind of PeriodicBoundaryIndex or via a more generic index with a boundary='periodic' option.

2 Likes
11

Hi! This was a great thread on handling circular lon values, and I wonder if anyone had experience on doing this with curvilinear grids? I’m working with the MPI-ESM1-2-LR CMIP6 model and I’m trying to subset patches of the N Pacific ocean that are located at the edges of the lat/lon grid of this model:

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(cat_url)
cat = col.search( table_id='Omon',experiment_id='historical',
             variable_id=['thetao'], grid_label='gn', member_id='r1i1p1f1',source_id='MPI-ESM1-2-LR') 
cmip6_compiled = cat.to_dataset_dict(
        zarr_kwargs={'consolidated':True, 'decode_times': True, 'use_cftime': True},
        aggregate=False)    
ds = xr.merge([v for k,v in cmip6_compiled.items()])
ds_subset = ds.where((ds.latitude > 0) & (ds.latitude < 60), drop=True).squeeze().load()
ds_subset['thetao'].isel(lev=0,time=0).plot(yincrease=False)

0 to 60 N

da = ds_subset.where((ds_subset.longitude > 125) & (ds_subset.longitude < 250), drop=True).squeeze().load()
da['thetao'].isel(lev=0,time=0).plot(yincrease=False)

N Pacific

I’ve used the xr.roll code mentioned above and it executes, however I get the ff error when I try to plot the resulting dataset:

ValueError: The input coordinate is not sorted in increasing order along axis 0. This can lead to unexpected results. Consider calling the `sortby` method on the input DataArray. To plot data with categorical axes, consider using the `heatmap` function from the `seaborn` statistical plotting library.

I’ve tried looking over the new functionality for handling these types of cases for xarray but I’m confused as to how to apply it to my dataset. I would appreciate any suggestions, thank you so much!