QGIS and grid_mapping attribute

[guziy]

Hi:

I noticed that if a converted field has grid_mapping attribute, qgis does not recognise the CRS, but if from the same file and variable I remove the grid_mapping attribute, then the CRS/projection is recognized by qgis... It is important for Geomet to accept the netcdf files that the CRS is correctly recognized by the tools like QGIS or gdal.

Here are some metadata of my example file:

# CRS not recognized for ETAS

$ ncdump -h 2019090606_.nc
netcdf \2019090606_ {                                                                                                                                                                                              
dimensions:                                                                                                                                                                                                        
        lat = 541 ;                                                                                                                                                                                                
        lon = 841 ;                                                                                                                                                                                                
        time = UNLIMITED ; // (6 currently)                                                                                                                                                                        
        sfclevel = 1 ;                                                                                                                                                                                             
variables:                                                                                                                                                                                                         
        float lat(lat) ;                                                                                                                                                                                           
                lat:long_name = "latitude" ;                                                                                                                                                                       
                lat:standard_name = "latitude" ;
                lat:units = "degrees_north" ;
        float lon(lon) ;
                lon:long_name = "longitude" ;
                lon:standard_name = "longitude" ;
                lon:units = "degrees_east" ;
        char crs_latlon ;
                crs_latlon:grid_mapping_name = "latitude_longitude" ;
                crs_latlon:earth_radius = 6371229. ;
        float ETAS(time, sfclevel, lat, lon) ;
                ETAS:_FillValue = 1.e+30f ;
                ETAS:coordinates = "leadtime reftime" ;
                ETAS:grid_mapping = "crs_latlon" ;
        double time(time) ;
                time:standard_name = "time" ;
                time:long_name = "Validity time" ;
                time:axis = "T" ;
                time:units = "hours since 2019-09-06 04:00:00" ;
                time:calendar = "gregorian" ;
        float sfclevel(sfclevel) ;
                sfclevel:axis = "Z" ;
                sfclevel:standard_name = "model_level_number" ;
                sfclevel:positive = "down" ;
        float leadtime(time) ;
                leadtime:standard_name = "forecast_period" ;
                leadtime:long_name = "Lead time (since forecast_reference_time)" ;
                leadtime:units = "hours" ;
        double reftime ;
                reftime:standard_name = "forecast_reference_time" ;
                reftime:units = "hours since 2019-09-06 03:00:00" ;
                reftime:calendar = "gregorian" ;

// global attributes:
                :Conventions = "CF-1.6" ;
                :history = "2021-09-16 21:04:32: /fs/ssm/eccc/crd/ccmr/EC-CAS/master/fstd2nc_0.20210305.1_all/bin/fstd2nc 2019090606_ 2019090606_.nc" ;
}


# CRS recognized

netcdf \2019090606_modified3 {
dimensions:
        lat = 541 ;
        lon = 841 ;
        time = UNLIMITED ; // (6 currently)
        sfclevel = 1 ;
variables:
        float lat(lat) ;
                lat:long_name = "latitude" ;
                lat:standard_name = "latitude" ;
                lat:units = "degrees_north" ;
        float lon(lon) ;
                lon:long_name = "longitude" ;
                lon:standard_name = "longitude" ;
                lon:units = "degrees_east" ;
        char crs_latlon ;
                crs_latlon:grid_mapping_name = "latitude_longitude" ;
                crs_latlon:earth_radius = 6371229. ;
        float ETAS(time, sfclevel, lat, lon) ;
                ETAS:_FillValue = 1.e+30f ;
                ETAS:coordinates = "leadtime reftime" ;
        double time(time) ;
                time:standard_name = "time" ;
                time:long_name = "Validity time" ;
                time:axis = "T" ;
                time:units = "hours since 2019-09-06 04:00:00" ;
                time:calendar = "gregorian" ;
        float sfclevel(sfclevel) ;
                sfclevel:axis = "Z" ;
                sfclevel:standard_name = "model_level_number" ;
                sfclevel:positive = "down" ;
        float leadtime(time) ;
                leadtime:standard_name = "forecast_period" ;
                leadtime:long_name = "Lead time (since forecast_reference_time)" ;
                leadtime:units = "hours" ;
        double reftime ;
                reftime:standard_name = "forecast_reference_time" ;
                reftime:units = "hours since 2019-09-06 03:00:00" ;
                reftime:calendar = "gregorian" ;

// global attributes:
                :Conventions = "CF-1.6" ;
                :history = "Thu Sep 16 21:24:58 2021: ncatted -a grid_mapping,ETAS,d,, 2019090606_.nc 2019090606_modified3.nc\n2021-09-16 21:04:32: /fs/ssm/eccc/crd/ccmr/EC-CAS/master/fstd2nc_0.20210305.1_all/bin/fstd2nc 2019090606_ 2019090606_.nc" ;
                :NCO = "4.7.2" ;

BTW is there a simple way to remove the sfclayer dimension during conversion?

[deacud]

Hi,

GDAL may still have problems properly interpreting cf-compliant metadata. As for Geomet, I remember that netcdf files generated by GIOPS-F could not be directly ingested into it for all sorts of reasons in the past. I am no longer up to date with such issues but will have to deliver netcdf files that Geomet is happy with in the coming weeks.

https://wiki.cmc.ec.gc.ca/wiki/GeoMet_2/Data/netCDF

https://gdal.org/drivers/raster/netcdf.html

"Georeference

There is no universal way of storing georeferencing in NetCDF files. The driver first tries to follow the CF-1 Convention from UNIDATA looking for the Metadata named “grid_mapping”. If “grid_mapping” is not present, the driver will try to find an lat/lon grid array to set geotransform array. The NetCDF driver verifies that the Lat/Lon array is equally space.

If those 2 methods fail, NetCDF driver will try to read the following metadata directly and set up georeferencing. ... "

[alexandreleroux]

Thanks @deacud . This page on the internal wiki is pretty old (2016!). GeoMet indeed leverages GDAL & PROJ to read input NetCDF files. The best way to test a NetCDF dataset is to provide a sample and we'll be glad to validate how it behaves. Thanks!

[neishm]

I'm not very familiar with the grid_mapping defaults (@deacud helped with those conventions). If I understand this thread correctly, grid_mapping is causing some issues with QGIS. However, it works ok with the GDAL & PROJ packages, so it should be ok for the intended purpose (using with GeoMet)? Or should there be a command-line option to remove it?

To remove the "sfclevel" dimension, you could pass the argument --exclude=sfclevel

[deacud]

For Sacha's case, regular lat-lon, 'grid_mapping' is necessary only if there is a need to specify the radius of the Earth. Otherwise, 'grid_mapping' could be removed.

"A grid mapping variable may be referenced by a data variable in order to explicitly declare the coordinate reference system (CRS) used for the horizontal spatial coordinate values. For example, if the horizontal spatial coordinates are latitude and longitude, the grid mapping variable can be used to declare the figure of the earth (WGS84 ellipsoid, sphere, etc.) they are based on."

"The grid_mapping variable may identify datums (such as the reference ellipsoid, the geoid or the prime meridian) for horizontal or vertical coordinates. Therefore a grid mapping variable may be needed when the coordinate variables for a horizontal grid are longitude and latitude. The grid_mapping_name of latitude_longitude should be used in this case."

https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html#grid-mappings-and-projections

[guziy]

Thanks, Daniel and Michael:

I am not sure how easy it is to add the command line option to remove grid_mapping, but I can still work around it by using ncatted from NCO and remove it after conversion. I thought maybe we (fstd2nc) needed to add more info to the longitude_latitude variable to make it work, that is mainly why this issue was created. If there is no other way but to remove the attribute, I can close the issue.

It would not be a bad thing if qgis could recognize the projection, if possible)

Cheers

[deacud]

"GDAL is a great library and a lot of work has gone into it, but its netcdf support has always been sketchy . When I first was directed to it years ago, it could only do 2-D files, and would flip the data, even when the metadata clearly said the axes went in the other direction (it just ignored the metadata attributes). That problem lasted for a long time (for all I know it still does this). GDAL has had problems with greater than 3-D files, forecast files, DSG files, files that are part of the NCEI examples for sending in data, some issues with time, and some of the newer features in netcdf4 files." (December 2019)

cf-convention/cf-conventions#222

[neishm]

It turns out to be trivial to remove the grid_mapping attribute (just need a small tweak to the --exclude code). However, as a side-effect it also removes the crs_latlon variable from the netcdf file. If that's acceptable I can commit the relevant changes.

[guziy]

It works when I remove crs_latlon, but I don't feel comfortable removing this information, especially if there are users out there relying on it (the variable and the attribute). I will try to see if I can get more info from qgis project, maybe the issue is there...

[neishm]

I tried replicating on my end, but not getting any error message with my test data. From ppp4, what I ran was:

fstd2nc /space/hall4/sitestore/eccc/crd/ccmr/min000/fstd2nc-testfiles/2017072800_ -f test.nc --vars=BTMY

This gives me a variable with the attribute grid_mapping = "cra_latlon". Then if I run

qgis test.nc

It brings up the data without any obvious problems.

@guziy, do you have a minimal working example of how to trigger the error?

edit: or is the problem that it doesn't pop up with a CRS dialogue for the file?
edit edit: for my file, there was no CRS dialogue when I opened it, but if I check the properties of the layer it does show CRS information

[guziy]

QGIS does open the file and shows the colors for data except in the lower right corner it says unknown CRS, so if you decide to use it along with shapefiles or other data, that might be problematic.

My example file is here:

/home/olh001/data/ppp3/sample_outputs/rdsps_v1.7.0/pseudo-analysis/2019090606_.nc

It is netcdf4 so you might want to use a newer qgis, from here for example:

/home/olh001/data/ords/miniconda3-backend/envs/qgis/bin/qgis

Let me know if you need more info.

[neishm]

Thanks. I was able to reproduce the same "Unknown CRS" message.

Clicking on "Unknown CRS" pulls up a mini map with the correct region selected. Also, overlaying a shapefile seems to give the right results. Maybe "Unknown" is a misnomer, and means more like "user-specified projection"?

For the default version of QGIS on ppp4 (2.18) says it has CRS of "USER:100000", and seems to work with a shapefile as well. I'm not familiar with GIS software... are there some other operations where this file fails to work?

[guziy]

No, I don't think other operations will fail if the ones you mention work.. I was a bit surprised by the fact that if I add the grid_mapping attribute it becomes unknown, but if it is a misnomer, then I guess this issue is a false alarm...

[neishm]

Out of curiosity, did the files end up working properly with QGIS in practical usage (beyond the quick checks that were done earlier)? Or is there further work needed to resolve this issue?

[guziy]

Thanks Mike, actually the files were accepted by GeoMet team, so all should be fine. I will close the issue, please, feel free to re-open if needed.