Shoud the nc file add CRS according to cf-conventions?

[singledoggy]

Recently I've learned the cf-conventions that defines a way to store crs in .nc file and pyproj gives a method to get crs from and to nc.
So is it possible to add crs when converting .nc? Otherwise we still need the original .ctl files to store the projection information. It's weird.
Thanks again for your work. It helps a lot in my work.

[miniufo]

That is interesting. Could you provide a more concrete example on how to do this? I just go through the doc you send, and the simplest way is to convert the crs in .ctl to a string of comments that will be put into the .nc? Then it is easy to convert it back using pyproj. Am I right?

[singledoggy]

I'm not a specialist of that too, so I guess you're right. And I don't know the details of this conventions and if it will become popular soon. As I know, many .ncfiles still don't save the projection information now. I just want to save compelete informations in xr.DataArray or the .nc files.

Xarray mentioned CF data just here, and there are few informations.

But this snippet returens True so at least this string can be safely converted back?

from pyproj import CRS
from xgrads import CtlDescriptor, get_data_projection
ctl = CtlDescriptor(file=ctl_path)
crs = get_data_projection(ctl)
cf = crs.to_cf()
crs_new = crs.from_cf(cf)
crs_new == crs

It seems that .nc files use an attribute named grid_mapping to store projection information.

[miniufo]

crs = get_data_projection(ctl)
cf = crs.to_cf()

get_data_projection does not return a CRS but a <cartopy.crs.LambertConformal at 0x11920270>. So I cannot find to_cf(). Have you tried this snippet?

[singledoggy]

I test that code in an environment contains cartopy==0.20.1, but in another version cartopy==0.17 or cartopy==0.19 it returns the same as your comment.

It seems that cartopy 0.20+ is required, and

cartopy.crs.CRS inherits from pyproj.crs.CRS, so it should behave like a pyproj.crs.CRS.

It may be a pretty new feature? But I know it sometimes cause compatibility issues.

[miniufo]

I guess this is really new feature. The doc shows that api's name is cs_to_cf() at here.

I currently cannot make any test because I don't have the newest environment. From the codes you showed it is pretty simple to do this. If cf is a string, adding it to the dataset before writting it to a disk file would be OK.

It seems that everything is changing to take into account the CF convention, although I knew it many years ago.

[miniufo]

When I tried to solve #37, I found a workaround for this. I use cartopy 0.18.0 and the returned crs from get_data_projection is not a subclass of pyproj.crs.CRS. But we can create one by:

from pyproj import CRS

crs = get_data_projection(ctl, Rearth=6370000)
proj_crs = CRS(crs.proj4_init)
cf = proj_crs.to_cf()

I get cf as:

It is a little too long and not a string. But crs.proj4_init is a short string as:

So if you want to put cf into nc file, you can put crs.proj4_init in a grid_mapping comment, and recover a CRS from the above code. I guess that cartopy is doing a significant update (I re-install anaconda and all python package and still get cartopy 0.18.0). We shall wait to see if the API's name is eventually fixed (e.g., the link you provided on 24 Mar at the very beginning is now broken...).

[singledoggy]

Yes, I've found it's an unsolved issue of Cartopy for a long time. And Xarray don't read the crs data by default. I've tried rioxarray to get the crs, but there are still some bugs.
So maybe this is their work to get a API to convert it.

[Mo-Dabao]

CF-Conventions are annoying but quite useful. I figured out the way to make panoply recgnize the crs in nc file:

https://bitegrads.readthedocs.io/en/latest/Examples/02_gridded_binary_with_pdef.html

The key here is to add an earth_radius in crs's attributes, otherwise panoply won't recgnize it:

https://github.com/Mo-Dabao/BiteGrADS/blob/ac279c3c4e0f6d2107b017a17ca11ff73ad266fe/src/bite_grads/grads.py#L157

[miniufo]

Hi @Mo-Dabao, that's very useful. Thank you for sharing this point. Recently, I've added a function to parse the global comments in ctl files output by WRF (see here).

I may add a util function in utils.py as (pseudo codes):

def to_CFnetcdf(dset, comments, Rearth, path):
    # add_comments_to_dset(dset, comments)
    # add_Rearth_to_dset(dset, Rearth)
    # dset.to_netcdf(path)

Do you think this is OK? Can you please try if the output NC file can be recognized by panoply, as I don't have panoply?

[Mo-Dabao]

I may add a util function in utils.py as (pseudo codes):

def to_CFnetcdf(dset, comments, Rearth, path):
    # add_comments_to_dset(dset, comments)
    # add_Rearth_to_dset(dset, Rearth)
    # dset.to_netcdf(path)

Do you think this is OK?

It's totally up to you. But the earth_radius I mentioned before may only apply to those data based on "spherical" Earth, for example some WRF outputs. So it might take more design to be compatible to all kinds of Earth.

Can you please try if the output NC file can be recognized by panoply, as I don't have panoply?

panoply can be download free from NASA's web site here. And it runs without installing, as long as you have JAVA 9 installed first.