2

I have two datasets:

Climate:

  • avg one month temperature
  • NetCDF File (CF Convention 1.6)
  • grid_mapping_name: rotated_latitude_longitude (with rlat and rlong coordinates defining the grid)
  • curvilinear grid

Harvest Area Fraction:

  • Wheat
  • Geotiff File
  • EPSG: 4326 (lon and lat coordinates defining the grid)
  • rectilinear grid

I managed to plot both datasets (see picture below).

I want to extract the climate data values where the wheat grows. However, the two datasets have different coordinate systems (rotated lon/lat grid vs. regular lon/lat grid), so they do not align.

Suggestions?

enter image description here

EDIT: If I call climatedata.rio.crs I get the following output:

GEOGCRS["undefined",BASEGEOGCRS["undefined",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION["Pole rotation (netCDF CF convention)",METHOD["Pole rotation (netCDF CF convention)"],PARAMETER["Grid north pole latitude (netCDF CF convention)",39.25,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["Grid north pole longitude (netCDF CF convention)",198,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["North pole grid longitude (netCDF CF convention)",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]

EDIT:

Step 1: Transform climate data to regular spaced lon-lat grid.

climatedata =  xa.open_dataset(r'filepath.nc', decode_cf = True, decode_coords = "all")
df = climatedata.squeeze().to_dataframe().reset_index()
geometry = gpd.points_from_xy(df.lon, df.lat)
gdf = gpd.GeoDataFrame(df, crs=climatedata.rio.crs, geometry=geometry)
geo_grid = make_geocube(vector_data=gdf, resolution=(-0.1, 0.1), rasterize_function=rasterize_points_griddata,)
geo_grid = geo_grid.tas_moy[:,:]

enter image description here

geo_grid.rio.crs is the same as mentioned above. harvestdata.rio.crs is the following:

GEOGCS["unknown",DATUM["unknown",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["unknown",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]

Step 2: Reproject climate data to harvest data:

harvestdata = xa.open_dataset(r'filepath.tif')
climatedata_matched = geo_grid.rio.reproject_match(harvestdata)

Problem: climatedata_matched is an array only with NaN values. enter image description here

Improve this question
3
  • What coordinate system is your climate data in?
    – Binx
    Commented Mar 31, 2022 at 14:10
  • The climate data has no 'crs' attribute. If I call 'climatedata.rio.crs', the output is posted above. I already tried to important an CRS object from the CF convention, but it returned no epsg code to me link. Is there another why to get the coordinate system?
    – ProgMau
    Commented Mar 31, 2022 at 14:57
  • Is it possible to share your files?
    – Binx
    Commented Mar 31, 2022 at 15:00

1 Answer 1

Reset to default
2

Step 1: Determine the CRS & coordinates for each one

For the climate grid with 2D coordinates, this may be a helpful reference: https://github.com/corteva/rioxarray/issues/209

Step 2: Re-project to the same grid

https://corteva.github.io/rioxarray/stable/examples/reproject_match.html

I recommend re-projecting the climate grid to the harvest grid.

Improve this answer
11
  • Very helpful post #209! Problem: The new matched climate grid only contains NaN values. Any suggestions? I will add my code to the question above.
    – ProgMau
    Commented Apr 1, 2022 at 12:15
  • You could try passing in the harvest grid into like= when making the geocube.
    – snowman2
    Commented Apr 2, 2022 at 3:43
  • Passing like=harvest data works. But the new plot has "stripes" meaning that the values for one x are the same along the y axis. Is this due to the rasterize function? That's why I tried rasterize_function=partial(rasterize_points_griddata, method="cubic") to keep the high resolution from the original data, but again I only get NaN values.
    – ProgMau
    Commented Apr 2, 2022 at 22:51
  • Afer analysing the data, my assumption is that reproject_match cannot read geo_grid.rio.crs (see output above). If I write a crs attribute manually to the geo_grid (geo_grid.rio.write_crs("epsg:4326", inplace=True)) then the reproject_match seems to work, because the array has values (not only NaN). Now my question is, how do I get the right epsg code from the climate data (netcdf)? Epsg:4326 was just a random assumption to try it out.
    – ProgMau
    Commented Apr 3, 2022 at 16:52
  • I am still not sure how to convert the wkt from climatedata.rio.crs (see above) to an epsg code. climatedata.rio.crs.to_epsg() returns none. Any suggestions on how to get the epsg code from the netcdf climate data?
    – ProgMau
    Commented May 15, 2022 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.