Keywords: ERA5, microsoft planetary computer, antimeridian, xarray, sel, sortby, stac search API
This post is motivated by playing around with the ERA5 data as hosted on Microsoft’s Planetary Computer.
I note that longitude in this data is over the range [0, 360), with 0 being the Greenwich Prime Meridian.
I also note that the STAC API for the bbox query parameter states
The coordinate reference system of the values is WGS 84 longitude/latitude (http://www.opengis.net/def/crs/OGC/1.3/CRS84)
Now, as I understand it, WGS84 is not a CRS per se, but describes the shape of the geoid upon which a particular CRS is based? Anyhoo, the STAC API documentation goes on to say
For WGS 84 longitude/latitude the values are in most cases the sequence
of minimum longitude, minimum latitude, maximum longitude and maximum
latitude. However, in cases where the box spans the antimeridian the
first value (west-most box edge) is larger than the third value
(east-most box edge).
So, here we’re effectively dealing with a dataset that plonks a/the antimeridian smack through the UK, which is kinda inconvenient for those of us who live here as opposed to vast swaths of the Pacific Ocean.
All of which is preamble, and arguably not of major relevance to the ERA5 data that is organised with one file for the whole Earth’s surface at a moment in time (as opposed to needing to locate which tiles/granules you need). Indeed, using bbox in the stac search here is pointless. This was my first learning. (Note, it wasn’t my first learning in time, but it’s worth sorting linearly with respect to the narrative here).
So I search for available data for my year of interest and don’t bother with the bbox parameter in the catalog search:
search = catalog.search(collections=['era5-pds'], datetime='2020', query={"era5:kind": {"eq": "fc"}})
datasets = [
xr.open_dataset(asset.href, **asset.extra_fields["xarray:open_kwargs"])
for item in search.items() for asset in item.assets.values()
]
ds = xr.combine_by_coords(datasets, join="exact")
ds
and get
Now, this appeals to me because I’ve got an xarray Dataset comprising dask.arrays of various interesting variables. But note the
lon goes from 0.0 to 359.8. I want to visualize Great Britain (sorry, N. Ireland) and xarray’s .sel() selector won’t give me what I want here. My next step is to transform those western longitudes:
ds['lon'] = ds.lon.where(ds.lon < 180, ds.lon - 360)
ds = ds.sortby('lon')
ds.lon
Okay, so now I have longitude running [-180, 180) with UK running from negative (west) to positive (east) and I can access using
.sel() and visualize:
uk_era5 = ds.sel(lat=slice(60, 50), lon=slice(-6, 2))
ax = uk_era5.air_temperature_at_2_metres_1hour_Maximum.groupby('time.month').min().sel(month=1).plot()
eng.plot(ax=ax.axes, facecolor='None')
which looks like I have the land in the right place.
Rounding off this post with my learnings/questions:
- Obviously, when accessing someone’s STAC API, you can’t change stuff server side, so when querying over the antimeridian, do you just need to break your query into two?
- And you need to know the coordinate system of the dimensions, because relying on the data complying with the STAC API and being ‘WGS84’ isn’t enough, because, like here, you can’t assume longitude runs [-180, 180)?
- Is it a good approach to transform the dimension coordinates (here to have negative values for the western hemisphere) and sort to ensure the coordinates monotonically increase?
Gently and politely pondering in the general direction of @TomAugspurger to ask whether they considered doing a coordinate (longitude) transformation of the data when ingesting? Because it strikes me as slightly odd that the STAC API seems silent on the actual CRS used. What’s to stop someone going “Well our data is referenced to the longitude of Time Square, New York, but it’s still on WGS84”? Is this a pointless question because in such cases where this matters to what you get back from the STAC API search (i.e. when the Earth’s surface is tessellated by tiles/granules) you’ll always be dealing with a longitude that runs [-180, 180)? Or do you just simply need to know the CRS of the data you’re querying in order to adjust or transform the bounds of your query parameters before making the call to search the catalog?