Thank you Even and Richard, Since Even opened a pull request on gdal's github ( https://github.com/OSGeo/gdal/pull/8407), I will continue the discussion over there.
Thomas lør. 16. sep. 2023 kl. 17:10 skrev Even Rouault <[email protected] >: > 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. > >
_______________________________________________ QGIS-User mailing list [email protected] List info: https://lists.osgeo.org/mailman/listinfo/qgis-user Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user
