<div dir="ltr">Thank you Even and Richard,<div><br></div><div>Since Even opened a pull request on gdal's github (<a href="https://github.com/OSGeo/gdal/pull/8407">https://github.com/OSGeo/gdal/pull/8407</a>), I will continue the discussion over there.</div><div><br></div><div>Thomas</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">lør. 16. sep. 2023 kl. 17:10 skrev Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Thomas,<br>
><br>
> I am not from the GIS world and not even a QGIS user, so bear with me <br>
> if my answers do not use the expected vocabulary.<br>
It takes years to get fully familar with GIS oddities, and netCDF is an <br>
area where even a life long of experience will not be enough to fully <br>
overcome the "creativity" of the netCDF community...<br>
><br>
> From my point of view, my netCDF files describe their geolocation in a <br>
> CF-compliant manner, with a grid_mapping CRS variable, and the x and y <br>
> projection coordinate variables. But QGIS does not recognize the CRS, <br>
> seemingly because I have my x and y as km (and I do specify <br>
> :units="km", I am not hiding this).<br>
><br>
> Since I am also a programmer, I tried to find where in the QGIS code <br>
> the netCDF/CF geolocation is read and decoded, to see if there was a <br>
> strict test on :units="m". But I failed to find this in the code.<br>
<br>
This is actually not really done in QGIS, but in PROJ (more exactly <br>
QGIS's QgsProjUtils::identifyCrs() delegates to PROJ's proj_identify() <br>
in <br>
<a href="https://github.com/qgis/QGIS/blob/8723b82ec19a3bec7f39b46128c798c7a2ee4230/src/core/proj/qgsprojutils.cpp#L271" rel="noreferrer" target="_blank">https://github.com/qgis/QGIS/blob/8723b82ec19a3bec7f39b46128c798c7a2ee4230/src/core/proj/qgsprojutils.cpp#L271</a> <br>
to try to correlate the CRS definition from the netCDF file to one in <br>
the database).<br>
<br>
As you likely open it as a raster (and not a mesh as Richard <br>
mentionned), the GDAL netCDF driver is also involved since it is it that <br>
reconstructs a CRS object from the attributes of the netCDF CF conventions.<br>
<br>
 From what I can tell, both PROJ and GDAL do the job "as expected". It <br>
is just that the way this CRS is encoded is not going to find any match <br>
with a known CRS.  To get what you want, the GDAL netCDF driver should <br>
both modify the CRS definition to change the km unit to m, and alter the <br>
geotransform matrix (which gives the coordinate of the top left corner <br>
and pixel size) to multiply their value by 1000. This could be done, but <br>
should it be done...? I'm not so sure. The issue is more than the <br>
"philosophy" of netCDF data producers is sometimes at odds with the <br>
usual CRS definitions.<br>
<br>
That said, for that very particular case, there's a proj4_string <br>
attribute in Lambert_Azimuthal_Equal_Area variable, with value <br>
"+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m <br>
+no_defs +type=crs", which indicates the intent to have a metre unit. <br>
There's also a global attribute geospatial_bounds_crs = "EPSG:6931" <br>
which further correlate this. Hence this enhancement to the GDAL netCDF <br>
driver to renormalize to metre: <a href="https://github.com/OSGeo/gdal/pull/8407" rel="noreferrer" target="_blank">https://github.com/OSGeo/gdal/pull/8407</a><br>
<br>
Even<br>
<br>
-- <br>
<a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
My software is free, but my time generally not.<br>
<br>
</blockquote></div>