Convert 2D gridded data from latlong projection in degrees to LAEA in meters

[Marston]

Hi,
I would like to take a 2D gridded dataarray with coordinates in degrees cover the globe and convert it to a LAEA projection in meters. I would also like to visualize the data both form in Xarray. But I'm not finding any good Xarray examples of this. I was looking at rasterio but that would mean that I would need to convert my gridded data to a raster image. Is there a better way to do this type of transformation in Xarray?

[dcherian]

at rasterio but that would mean that I would need to convert my gridded data to a raster image.

Do you have a vector dataset? "2D gridded datarray" sounds pretty raster-ish to me. It would help if you detailed your concerns with rasterio.

[Marston]

Here's an example of the vector dataarray:

import xarray as xr

# Load the DataArray from the NetCDF file
da = xr.open_dataset('gfs_orography.nc', engine='netcdf4').to_array()

# Define the projection parameters for the laea projection
projection_params = {
    'proj': 'laea',
    'lat_0': 0,
    'lon_0': -135,
    'a': 6370997.0,  # Semi-major axis of the ellipsoid
    'b': 6370997.0,  # Semi-minor axis of the ellipsoid
    'units': 'm',
    'no_defs': True
}

The idea here is to convert da from a EPSG:4326 to a LAEA defined globally so that I can then take a slice out the region as defined by projection_params.

[scottyhq]

@Marston if you are working with projections (e.g. anything PROJ can handle https://proj.org/en/9.3/operations/projections/laea.html) look towards rioxarray. Here is a similar discussion that might help, corteva/rioxarray#583.

[Marston]

I'm looking at rioxarray but I cannot find a viable solution. Here's my code:

import xarray
import rioxarray # for 'rio' accessor
import json as js
import pyproj
from pyproj import CRS

# Load region proj string from a JSON file
with open('domains.json', 'r') as json_file:
    region_config = js.load(json_file)

# Read and curate netcdf file
crs = CRS(init='EPSG:4326')
da = xarray.open_dataarray('gfs_orography.nc', engine='netcdf4', decode_coords='all')
# Rename
da = da.rename({'lon':'x', 'lat':'y'}).transpose('y', 'x')
da.rio.write_crs(crs.to_proj4(), inplace=True)


region = 'roi'
info = region_config[region]
# Extract parameters from the definition
x_len = info["xn"]
y_len = info["yn"]
projString = info['projection']
# Create a Pyproj CRS object from the proj string
crs = pyproj.CRS(projString)
# Get the EPSG code of the CRS
epsg_code = crs.to_epsg()
if epsg_code is None:
    print("EPSG Code:", epsg_code)
    print('Setting to: EPSG:9820')
    epsg_code = "EPSG:9820"
    # Calculate coordinates for the sector polygon
    x_min, x_max = -x_len / 2, x_len / 2
    y_min, y_max = -y_len / 2, y_len / 2
else:
    refloc = info['refoc']
    # Extract coordinates for the sector polygon
    x_min, x_max = refloc['xRef'], abs(refloc['xRef'])
    y_min, y_max = refloc['yRef'], abs(refloc['yRef'])
resolution = info['bin']

# Define the spatial extent
extent = (x_min, y_min, x_max, y_max)
# Define the resolution in meters (x, y)
resolution_x = -resolution
resolution_y = resolution

da_laea = da.rio.reproject(epsg_code)

gives the following error:

Exception has occurred: CRSError
The WKT could not be parsed. OGR Error code 6
rasterio._err.CPLE_BaseError: OGR Error code 6

During handling of the above exception, another exception occurred:

  File "/root/odc_test.py", line 51, in <module>
    da_laea = da.rio.reproject(epsg_code)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^
rasterio.errors.CRSError: The WKT could not be parsed. OGR Error code 6

Been searching for a clue as to why this is failing but coming up short. I'm using the latest versions of xarray and rioxarray in python 3.11.