Convert ERA5 data to WGS84

Question

I'm using the cdsapi python package to download ERA5 data in netcdf format. I'd like to plot it alongside data from a GPS track and so I'm trying to ensure that the coordinate reference systems (CRS) of both datasets are the same. Ideally, I'd like to convert the ERA5 data to WGS84.

When I open the ERA5 netcdf data with xarray, no CRS is recognized. I've found the spatial reference/earth model used in ERA5 is a sphere with radius 6367.47 km (assuming netcdf inherits from GRIB1 instead of GRIB2), but I can't find a name for this (incomplete?) CRS and don't understand how to use this information to convert the ERA5 data to WGS84.

The comment on this question suggests that I can just convert the longitude values of my ERA5 data to values between -180 and 180, set the CRS to WGS84, and be done with it. But this doesn't seem right, because the WGS84 uses an ellipsoid rather than a sphere to derive the coordinates of physical locations. How can I best convert the ERA5 data to WGS84? I'm interested in both the relevant theory and its practical application (preferably in python).

Here's a minimum working example of my current approach. First I download data using the following code.

import cdsapi
client = cdsapi.Client()
client.retrieve(
    'reanalysis-era5-single-levels', {
        'product_type': 'reanalysis',
        'variable':     '2m_temperature',
        'year':         '2017',
        'month':        '12',
        'day':          '17',
        'time':         '00:00',
        'area':         [90, 0, 0, 360],
        'grid':         [1.0, 1.0],
        'format':       'netcdf'
},
'era5_test.nc')

Then I load the data with xarray and convert it to a GeoDataFrame.

import xarray as xr
import geopandas as gpd
ds = xr.open_dataset('era5_test.nc')
df = ds.to_dataframe().reset_index()
gdf = gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df['longitude'],
                                                         df['latitude']))

Then I transform the longitude values and set the CRS to WGS84.

import numpy as np

def transform_era5_lon_to_wgs84_lon(lon):
    return np.mod(lon + 180, 360) - 360

gdf = gdf.apply(lambda col: col.apply(tranfrom_era5_lon_to_wgs84_lon)
                if col.name == 'longitude' else col)
gdf.crs = {'init': 'epsg:4326'}

Answer 1

I ended up reaching out to the Copernicus support team and it turns out that ERA5 data is already referenced in the horizontal (lat/lon) with respect to the WGS84 ellipse so there is no need to convert it. They pointed me to this link (link dead, a similar sentence is now here), which includes the following statement:

ECMWF data is referenced in the horizontal with respect to the WGS84 ellipse (which defines the major/minor axes) but in the vertical it is referenced to the Geoid (EGM96).

Answer 2

If you are using a NetCDF file, then your data is in EPSG:4326 geographic coordinates.

To rotate the coordinates: https://github.com/corteva/rioxarray/issues/58

However, if you load in the native grib format, then you would want to build a custom geographic CRS:

https://pyproj4.github.io/pyproj/stable/build_crs.html#geographic-crs

I would recommend using rioxarray to write the CRS to the file: https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html

And to reproject: https://corteva.github.io/rioxarray/stable/examples/reproject.html

Answer 3

Using Juanma Cintas' solution I found a simpler way of "flipping" the variable in the Xarray dataset was to use the .sortby(coordinate) method after converting the 0-360 longitude to standard -180 - 180.

lonaxis = wd['longitude'].values
newlonaxis = lonaxis[lonaxis >= 180] -= 360

wd['longitude'] = lonaxis

wd.sortby('longitude').sel(time='2020-01-01T15:00:00', method='nearest')['tc_ch4'].plot()

Which yields a coherent data array

Answer 4

I've found myself in a similar situation, where I had an ERA5 image centered on the pacific since its coordinates were ranging from 0 to 360. I've based myself on all your comments and links, but I didn't fully get it without flipping the raster values too, not only the longitude coordinates. I will copy here my working solution at the moment in case someone else has the exact same problem:

def fromERA5toWGS84(dataarray, output=None, xaxis = "longitude", name = None):
    import rioxarray
    import xarray as xr
    if type(dataarray) is str:
        if dataarray.split(".")[-1] == "nc":
            dataarray = xr.open_dataset(dataarray, decode_coords="all")
        else:
            dataarray = xr.open_rasterio(dataarray, decode_coords="all")

    else:
        pass
    if type(dataarray) is xr.core.dataset.Dataset:
        if name is None: raise ValueError("I need a variable name")
        dataarray = dataarray[name]

    longitud = dataarray[xaxis].values
    longitud[longitud >= 180] -= 360
    indexes = np.where(longitud >= 0)[0].tolist()
    outindexes = np.where(longitud < 0)[0].tolist()
    array = dataarray.values
    zeros = np.zeros(array.shape)
    zeros[zeros==0] = np.nan
    flipped1 = array[:,:,indexes]
    flipped2 = array[:,:,outindexes]
    zeros[:,:,indexes] = flipped2
    zeros[:,:,outindexes] = flipped1
    longitud.sort()
    dataarray[xaxis] = longitud
    dataarray.values = zeros
    dataarray = dataarray.rio.write_crs(4326)
    if output is not None:
        if output.split(".")[-1] == "nc":
            dataarray.to_netcdf(output)
            del dataarray
            return "Done!"
        elif output.split(".")[-1] == "tif":
            dataarray.rio.to_raster(output)
            del dataarray
            return "Done!"
        else:
            raise ValueError("format not implemented")
    return dataarray