[Qgis-user] Does QGIS require units "m" for the projection coordinate variables in a netCDF/CF file?

[Even Rouault]

Thomas,
>
> I am not from the GIS world and not even a QGIS user, so bear with me 
> if my answers do not use the expected vocabulary.
It takes years to get fully familar with GIS oddities, and netCDF is an 
area where even a life long of experience will not be enough to fully 
overcome the "creativity" of the netCDF community...
>
> From my point of view, my netCDF files describe their geolocation in a 
> CF-compliant manner, with a grid_mapping CRS variable, and the x and y 
> projection coordinate variables. But QGIS does not recognize the CRS, 
> seemingly because I have my x and y as km (and I do specify 
> :units="km", I am not hiding this).
>
> Since I am also a programmer, I tried to find where in the QGIS code 
> the netCDF/CF geolocation is read and decoded, to see if there was a 
> strict test on :units="m". But I failed to find this in the code.

This is actually not really done in QGIS, but in PROJ (more exactly 
QGIS's QgsProjUtils::identifyCrs() delegates to PROJ's proj_identify() 
in 
https://github.com/qgis/QGIS/blob/8723b82ec19a3bec7f39b46128c798c7a2ee4230/src/core/proj/qgsprojutils.cpp#L271 
to try to correlate the CRS definition from the netCDF file to one in 
the database).

As you likely open it as a raster (and not a mesh as Richard 
mentionned), the GDAL netCDF driver is also involved since it is it that 
reconstructs a CRS object from the attributes of the netCDF CF conventions.

 From what I can tell, both PROJ and GDAL do the job "as expected". It 
is just that the way this CRS is encoded is not going to find any match 
with a known CRS.  To get what you want, the GDAL netCDF driver should 
both modify the CRS definition to change the km unit to m, and alter the 
geotransform matrix (which gives the coordinate of the top left corner 
and pixel size) to multiply their value by 1000. This could be done, but 
should it be done...? I'm not so sure. The issue is more than the 
"philosophy" of netCDF data producers is sometimes at odds with the 
usual CRS definitions.

That said, for that very particular case, there's a proj4_string 
attribute in Lambert_Azimuthal_Equal_Area variable, with value 
"+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m 
+no_defs +type=crs", which indicates the intent to have a metre unit. 
There's also a global attribute geospatial_bounds_crs = "EPSG:6931" 
which further correlate this. Hence this enhancement to the GDAL netCDF 
driver to renormalize to metre: https://github.com/OSGeo/gdal/pull/8407

Even

-- 
http://www.spatialys.com
My software is free, but my time generally not.