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 ice_2021_07.nc")
print(ds)
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 ...
shapefile
gdf = gpd.read_file('polygon.shp').geometry.to_list()
print(gdf)
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.
ds.rio.write_crs("epsg:4326", inplace=True)
clipped = ds.rio.clip(gdf, crs="epsg:4326")
print(clipped)
The output looks like :
<xarray.Dataset>
Dimensions: (xc: 13, yc: 1)
Coordinates:
* 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