Xarray and irregular grid (4 x 5 °, poles are different size)

Science
1

Hello,

I’m trying to create an xarray dataset with an irregular grid. It’s a 4x5 ° lon by lat (72 x 46). Except the polar bin is 2 rather than 4 degrees.

I can make a GeoPandas grid that matches this with the following code.

import geopandas as gpd
import shapely

xres, yres = 5,4
lon_bnds = np.arange(-180,180+1,xres)
lat_bnds = np.append(np.append(-90, np.arange(-88,88+1,yres)),90)
polygons = []
for i in range(lon_bnds.size-1):
    lon0,lon1 = lon_bnds[[i,i+1]]
    for j in range(lat_bnds.size-1):
        lat0,lat1 = lat_bnds[[j,j+1]]
        polygons.append(
            shapely.Polygon([(lon0,lat0), (lon1,lat0), (lon1,lat1), (lon0,lat1)])
        )

grid = gpd.GeoDataFrame({'geometry': polygons}, crs='epsg:4326')
grid['id'] = grid.index

sel = grid.intersects(shapely.Polygon([(-180,-90), (-100,-90), (-100,-72), (-180,-72)]))
_ = grid[sel].plot(color='none', edgecolor='gray')

And the graphic shows the S. Pole cell is 2 degrees (-88 to -90)

I’m not sure how to make the xarray equivalent of this. I’ve tried the follow. Plotting the bottom left corner, the bottom row starts at -87.5, not -88.

import numpy as np
import xarray as xr
import pandas as pd
# import rioxarray as rio

xres, yres = 5,4
lon_bnds = np.arange(-180,180+1,xres)
lat_bnds = np.append(np.append(-90, np.arange(-88,88+1,yres)),90)

lats = np.append(np.append(-89, np.linspace(-86,86,44)), 89)
lons = np.linspace(-177.5,177.5,72)

data = np.arange(lats.size*lons.size).reshape(lats.size,lons.size)%5

foo = xr.DataArray(
    data = data,
    dims = ['lat','lon'],
    coords = {'lat': lats, 'lon': lons}
)
ds = xr.Dataset({'foo':foo})
ds = ds.rio.write_crs('epsg:4326')
ds = ds.rio.set_spatial_dims('lon','lat')

ds['lat_bnds'] = (('lat_bnds'), lat_bnds)
ds['lon_bnds'] = (('lon_bnds'), lon_bnds)

ds['lat'].attrs['bounds'] = 'lat_bnds'
ds['lat'].attrs['long_name'] = 'latitude'

print(ds)
_ = ds['foo'].sel({'lat':slice(-90,-80), 'lon':slice(-180,-160)}).plot()

output:

    <xarray.Dataset>
    Dimensions:      (lat: 46, lon: 72, lat_bnds: 47, lon_bnds: 73)
    Coordinates:
      * lat          (lat) float64 -89.0 -86.0 -82.0 -78.0 ... 78.0 82.0 86.0 89.0
      * lon          (lon) float64 -177.5 -172.5 -167.5 -162.5 ... 167.5 172.5 177.5
        spatial_ref  int64 0
      * lat_bnds     (lat_bnds) int64 -90 -88 -84 -80 -76 -72 ... 72 76 80 84 88 90
      * lon_bnds     (lon_bnds) int64 -180 -175 -170 -165 -160 ... 165 170 175 180
    Data variables:
        foo          (lat, lon) int64 0 1 2 3 4 0 1 2 3 4 0 ... 2 3 4 0 1 2 3 4 0 1

image

2

I was only allowed to embed one graphic, so here is the grid from the first bit of code, which does show -88 to -90 cells.

image

3

Welcome @mankoff! :wave:

You’re doing everything right. The problem is that Xarray doesn’t actually use the bounds variables when plotting. Instead, it just makes up cell “boundaries” using some ad-hoc logic. Unlike GeoPandas, Xarray doesn’t explicitly model the geometry of “grid cells.”

However, this can and should change in Xarray. See this issue from @dcherian for example:

The cf-xarray package has some deeper support for the bound variables. (The way you are encoding your bounds is part of CF conventions.) - Bounds Variables - cf_xarray documentation

Finally, the xvec packages aims to add geopandas-style geometry support to Xarray:

https://xvec.readthedocs.io/en/stable/

For now, to get an exact plot, you can’t use Xarray. You have to just use lower-level matplotlib code, explicitly specifying the coordinate variables for x and y, i.e. something like

plt.pcolormesh(lat_bnds, lon_bnds, data_to_plot)
3 Likes
4

Ah nice, I was wondering :ok_hand::grin:

1 Like
5

I think for now I’ll work on a regular grid from -88 to 88, and then deal with the poles separately at the end.

Thanks for the info and encouragement that I’m doing something right :).

1 Like