Clipping an xarray Dataset (CDS Climate Sea Ice) using a Shapefile with non-standard grid coordinates


I have a netCDF dataset with dimensions represented by ‘xc’ and ‘yc’, which are in meters, and the CRS of the dataset is Lambert Azimuthal Equal Area projection with the following parameters:

Grid Mapping Name: Lambert Azimuthal Equal Area
Longitude of Projection Origin: 0.0
Latitude of Projection Origin: -90.0
False Easting: 0.0
False Northing: 0.0
Semi-major Axis: 6378137.0
Inverse Flattening: 298.257223563
Proj4 String: +proj=laea +lon_0=0 +datum=WGS84 +ellps=WGS84 +lat_0=-90.0

I also have a shapefile that represents the area I want to crop from the dataset.

I’m trying to crop this dataset using the shapefile, but I’m facing some difficulties due to the non-standard grid coordinates and CRS.

Here’s what I’ve tried so far:

import geopandas as gpd
import rioxarray
import xarray as xr

# Open the netCDF dataset
ds = xr.open_dataset("./sea")

Here is the output :

> <xarray.Dataset> Dimensions:   (xc: 432, yc: 432) Coordinates:   * xc 
> (xc) float64 -5.388e+03 -5.362e+03 ... 5.362e+03 5.388e+03   * yc     
> (yc) float64 5.388e+03 5.362e+03 ... -5.362e+03 -5.388e+03
>     time      datetime64[ns] ...
>     lat       (yc, xc) float32 ...
>     lon       (yc, xc) float32 ... Data variables:
>     ice_conc  (yc, xc) float64 ...


gdf = gpd.read_file('polygon.shp').geometry.to_list()

Shapefile output in usual lat long format:

[<POLYGON ((-180 -64.4, -180 -64.4, -180 -64.5, -180 -64.5, -180 -64.6, -180 ...>]

The problem occurred when running the code rio.clip. It shows only one lat (should be more) and the cropped image is entirely null/empty."epsg:4326", inplace=True)
clipped =, crs="epsg:4326")

The output looks like :

Dimensions:                 (xc: 13, yc: 1)
  * xc                      (xc) float64 -137.5 -112.5 -87.5 ... 137.5 162.5
  * yc                      (yc) float64 -62.5
    time                    datetime64[ns] ...
    lat                     (yc, xc) float32 ...
    lon                     (yc, xc) float32 ...
    Lambert_Azimuthal_Grid  int64 0
Data variables:
    ice_conc                (yc, xc) float64 nan nan nan nan ... nan nan nan nan

Why this is happening? I think, there is a problem in non-standard grid coordinates.

Could someone guide me on how to properly crop the dataset using the shapefile, taking into account the non-standard grid coordinates and CRS?

Any help would be greatly appreciated!

Dataset and shapefile are here drive