Handling slicing with circular longitude coordinates in xarray

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

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

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

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

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

@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!

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

cc @benbovy

2 Likes

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)
2 Likes

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