reproject and clip netcdf file

[dhetting]

I have a netcdf file I'm trying to reproject and clip, but I'm having a hard time sorting out the georeferencing. Everything plots as expected in Panopoly.

When I load the file with ds = rioxarray.open_rasterio('daymet_v4_daily_na_prcp_2020.nc'), I get a warning: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs., which I take to mean the file is formatted improperly for rioxarray and is probably the source of my problems, but I'm not sure what to do to fix that.

My dataset:

I'm interested in the prcp values, which I see have the lambert_conformal_conic grid_mapping:

 'lambert_conformal_conic#false_northing': 0,
 'lambert_conformal_conic#grid_mapping_name': 'lambert_conformal_conic',
 'lambert_conformal_conic#inverse_flattening': 298.257223563,
 'lambert_conformal_conic#latitude_of_projection_origin': 42.5,
 'lambert_conformal_conic#longitude_of_central_meridian': -100,
 'lambert_conformal_conic#semi_major_axis': 6378137,
 'lambert_conformal_conic#standard_parallel': array([25., 60.]),
 'citation': 'Please see http://daymet.ornl.gov/ for current Daymet data citation information',
 'Conventions': 'CF-1.6',
 'references': 'Please see http://daymet.ornl.gov/ for current information on Daymet references',
 'source': 'Daymet Software Version 4.0',
 'start_year': 2020,
 'Version_data': 'Daymet Data Version 4.0',
 'Version_software': 'Daymet Software Version 4.0',
 'NETCDF_DIM_EXTRA': '{time}',
 'NETCDF_DIM_time_DEF': array([365.,   5.]),
 'NETCDF_DIM_time_VALUES': array([25567.5, ..., 25931.5]),
 'prcp#cell_methods': 'area: mean time: sum',
 'prcp#coordinates': 'lat lon',
 'prcp#grid_mapping': 'lambert_conformal_conic',
 'prcp#long_name': 'daily total precipitation',
 'prcp#missing_value': -9999,
 'prcp#units': 'mm/day',
 'prcp#_FillValue': -9999,
 'time#bounds': 'time_bnds',
 'time#calendar': 'standard',
 'time#long_name': '24-hour day based on local time',
 'time#standard_name': 'time',
 'time#units': 'days since 1950-01-01 00:00:00',
 'x#long_name': 'x coordinate of projection',
 'x#standard_name': 'projection_x_coordinate',
 'x#units': 'm',
 'y#long_name': 'y coordinate of projection',
 'y#standard_name': 'projection_y_coordinate',
 'y#units': 'm'}

When I try to reproject with ds[0].rio.reproject('epsg:4326'), it fails with MissingCRS: CRS not found. Please set the CRS with 'rio.write_crs()'. Data variable: lat.

So I go to write the CRS:

crs = CRS(ds[0].lambert_conformal_conic.crs_wkt)

ds[0].rio.write_crs(crs.to_string(), inplace=True)

but that fails with RioXarrayError: Multiple grid mappings exist..

How should I proceed in order to reproject my dataset? Thanks!

[snowman2]

The issue is likely to do with your lon/lat variables as they are 2D. I would select the variable(s) you want to use:

ds = rioxarray.open_rasterio('daymet_v4_daily_na_prcp_2020.nc', variable="prcp")

Also, if grid_mapping is already in your attrs/encoding then you don't need to write the CRS.

[snowman2]

You can also use xarray.open_dataset:

xds = xarray.open_dataset("daymet_v4_daily_na_prcp_2020.nc", decode_coords="all")

It seems to recognize the CRS just fine:

xds.rio.crs

CRS.from_wkt('PROJCS["undefined",GEOGCS["undefined",DATUM["undefined",SPHEROID["undefined",6378137,298.257223563]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",60],PARAMETER["latitude_of_origin",42.5],PARAMETER["central_meridian",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

[dhetting]

Great, that did it. Thanks for the quick reply!