netCDF: renormalize CRS and geotransform to metric, for EUMETSAT OSI …

[rouault]

…SAF products with a proj4_string, units=km and geospatial_bounds_crs

Cf https://lists.osgeo.org/pipermail/qgis-user/2023-September/053413.html

[coveralls]

coverage: 67.612% (+0.002%) from 67.61% when pulling 9e63a92 on rouault:netcdf_EUMETSAT_OSI_SAF into 43e51bc on OSGeo:master.

[TomLav]

Thanks @rouault for this quick commit.

I'd like to give some more context to avoid that you bring in code that are not strictly required.

First I want to mention that this is not only the OSI SAF netCDF files that can use "km" as a unit for the projection x/y coordinate variables. For good or bad, this is allowed in the netCDF/CF convention (see Example "5.7. Lambert conformal projection" in CF-1.10). By allowing and acting on :units = "km", gdal will support many CF-compliant datasets.

Second, from what I can see in the CF convention, there is nothing special about the ":units" attribute of the projection x/y coordinate variables (when not "degrees"). These are thus treated as any :units attribute, by the "udunit" package, which supports all unit suffixes from yotta: 1e24 to yocto: 1e-24 (See table 3.1. Supported Units in CF-1.10). Should gdal natively support files whose projection x/y coordinate variables are specified in Gm (megameters) or um (micrometers) I do not think so. But the "km" is part of a bigger picture.

IMO, only the the CF-compliant description of the CRS should be decoded by GDAL:

int Lambert_Azimuthal_Equal_Area ;
	Lambert_Azimuthal_Equal_Area:grid_mapping_name = "lambert_azimuthal_equal_area" ;
	Lambert_Azimuthal_Equal_Area:longitude_of_projection_origin = 0.f ;
	Lambert_Azimuthal_Equal_Area:latitude_of_projection_origin = 90.f ;
	Lambert_Azimuthal_Equal_Area:false_easting = 0. ;
	Lambert_Azimuthal_Equal_Area:false_northing = 0. ;
	Lambert_Azimuthal_Equal_Area:semi_major_axis = 6378137. ;
	Lambert_Azimuthal_Equal_Area:inverse_flattening = 298.257223563 ;

Note that in the above I removed the proj4_string attribute.

In the discussion at the osgeo users mailing list, you rightly pointed to at the "creativity" and "phylosophy" of the netCDF community. While I am certain that you and many others have had to add special treatments to support netCDF files (like this feature), I don't think it is a good idea to allow any strange idea to your tool.

For example, we have a "proj4_string" because of historical reasons (a tool we needed to use could not decode the CF description of the CRS natively) and to help human readers of our files who want an easy access to the proj string. We could easily remove this attribute from later version of our code because as far as we are concerned, the CF description of the CRS rules. We would not want to have to keep it so that gdal decodes our files.

We also have a ":geospatial_bounds_crs" global attribute because of the ACDD-1.3 convention (separate from CF but that often go in pair) that strongly recommend it.

You are free to bring your commit in full, but in my view the only required part is the handling of the "km" units for the projection x/y coordinate variables, so that gdal will support many CF-compliant datasets (including this OSI SAF data).

I do not think you should check on ":proj4_string" because this is not part of the CF standard. The case of ":geospatial_bounds_crs" is a bit tricky as it is Recommended by ACDD-1.3, but is not part of the CF convention. I do not think it is gdal's role to check if the CRS specified using the CF convention is the same as the CRS specified by the ACDD attribute.

I hope the above helps limiting the commit to the part that is required, and avoid bringing in code that is not strictly following the CF convention.

[rouault]

The thing is that GDAL already "properly" supported units=km on the X and Y axis, by incorporating it into the CRS definition, as this seems to be the one logical way of dealing with it. But when doing so, this prevents identification of the CRS with any EPSG registered CRS when units=km, but this is definetely the proper thing to do when units="us survey foot" (since a lot of US CRS are defined like that). So one could do the special processing for units=km to put the CRS axis unit to be metre, and multiply the axis values by 1000 in all cases (independently of the presence of proj4_string and geospatial_processing_bounds). I was just a bit worried about breaking backward compatibility where people would have processing chains expecting the CRS to be in km, but perhaps this is not such a bit concern and/or we could offer another open option (PRESERVE_AXIS_UNIT_IN_CRS=YES/NO).

@mdsumner any opinion ?

[rouault]

@mdsumner cf above (sorry if you get double pinged, but I had your github handle wrong the first time, and I'm not sure if one gets notification when fixing it in an existing comment)

[mdsumner]

all good appreciate the ping, and ongoing discussion 🙏 I think the preserve axis option is a good move

[rouault]

extra commit to make the change less specific to EUMETSAT OSI SAF and add a PRESERVE_AXIS_UNIT_IN_CRS=YES/NO open option, which defaults to NO

[TomLav]

Thank you for taking my input into account. I do not know enough about GDAL to assess if this commit triggers corner cases.

I can assist in testing this upgrade on more files once it gets into a gdal release (I wouldn't know how to do this on the github version). Does GDAL provide a (command line) tool to dump CRS + AXISCOORDs to screen given a filename?

[rouault]

I can assist in testing this upgrade on more files once it gets into a gdal release

In a few hours through GDAL master conda builds: https://gdal.org/download.html#gdal-master-conda-builds

Does GDAL provide a (command line) tool to dump CRS + AXISCOORDs to screen given a filename?

gdalinfo: https://gdal.org/programs/gdalinfo.html