Creating (x,y) coordinates in .nc from .grib2

Question

I have to use a .nc file as input in a model and such model requires projection_x_coordinate and projection_y_coordinate as variables but I don't have them. Is it possible to create them using pyproj, for instance? (file available here). My projparams are:

{'a': 6371229,
 'b': 6371229,
 'proj': 'stere',
 'lat_ts': 60.0,
 'lat_0': 90.0,
 'lon_0': 264.0}

My first idea was to use the information of the number of points (Nx in the file), the horizontal resolution (Dx) and the first points (e.g. latitudeOfFirstGridPointInDegrees) to create the desired variables but I don't know if it is the right thing to do.

snowman2 · Accepted Answer

You may be interested in rioxarray.
EDIT: with rioxarray>=0.0.31 and rasterio>=1.1.5 you should only need to do:
import rioxarray
rds = rioxarray.open_rasterio("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2")
rds.to_netcdf("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005_coord.nc")

Otherwise:
Due to https://github.com/mapbox/rasterio/issues/1248, you cannot use rioxarray.open_rasterio to generate CRS/coordinate info.
So, you will have to do it manually:
from rioxarray.rioxarray import affine_to_coords
import rasterio

with rasterio.Env(), rasterio.open("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2") as rrr:
    riocrs = rrr.crs
    coords = affine_to_coords(rrr.transform, rrr.width, rrr.height)

Then, you can use rioxarray to read in the data, add attrs, and write to netcdf:
from pyproj.crs import CRS
import rioxarray

rds = rioxarray.open_rasterio("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005.grib2", parse_coordinates=False)
rds = rds.assign_coords(coords)
rds.rio.write_crs(riocrs, inplace=True)
rds.spatial_ref.attrs.update(CRS.from_user_input(riocrs).to_cf())
rds.x.attrs["long_name"] = "x coordinate of projection"
rds.x.attrs["standard_name"] = "projection_x_coordinate"
rds.y.attrs["long_name"] = "y coordinate of projection"
rds.y.attrs["standard_name"] = "projection_y_coordinate"
rds.to_netcdf("CMC_caps_HGT_ISBL_0050_ps3km_2018090900_P005_coord.nc")

Creating (x,y) coordinates in .nc from .grib2

One Answer

Add your own answers!

Ask a Question