netCDF: use the SRS (its geographic part) if found in the file, inste…

…ad of the hardcoded WGS84 string, for the GEOLOCATION.SRS metadata item (fixes #9526)
OSGeo · Mar 22, 2024 · 5882838 · 5882838
1 parent d67906c
commit 5882838
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 8 deletions.
diff --git a/autotest/gdrivers/netcdf.py b/autotest/gdrivers/netcdf.py
@@ -1576,11 +1576,13 @@ def test_netcdf_42():
     gdaltest.netcdf_drv.CreateCopy("tmp/netcdf_42.nc", src_ds)
 
     ds = gdal.Open("tmp/netcdf_42.nc")
-    assert ds.GetMetadata("GEOLOCATION") == {
+    got_md = ds.GetMetadata("GEOLOCATION")
+    assert got_md["SRS"].startswith('GEOGCRS["WGS 84",')
+    del got_md["SRS"]
+    assert got_md == {
         "LINE_OFFSET": "0",
         "X_DATASET": 'NETCDF:"tmp/netcdf_42.nc":lon',
         "PIXEL_STEP": "1",
-        "SRS": 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]',
         "PIXEL_OFFSET": "0",
         "X_BAND": "1",
         "LINE_STEP": "1",
@@ -1718,11 +1720,13 @@ def test_netcdf_43():
     )
 
     ds = gdal.Open("tmp/netcdf_43.nc")
-    assert ds.GetMetadata("GEOLOCATION") == {
+    got_md = ds.GetMetadata("GEOLOCATION")
+    assert got_md["SRS"].startswith('GEOGCRS["NAD27",')
+    del got_md["SRS"]
+    assert got_md == {
         "LINE_OFFSET": "0",
         "X_DATASET": 'NETCDF:"tmp/netcdf_43.nc":lon',
         "PIXEL_STEP": "1",
-        "SRS": 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]',
         "PIXEL_OFFSET": "0",
         "X_BAND": "1",
         "LINE_STEP": "1",

diff --git a/frmts/netcdf/netcdfdataset.cpp b/frmts/netcdf/netcdfdataset.cpp
@@ -4266,9 +4266,25 @@ void netCDFDataset::SetProjectionFromVar(
     if (bReadSRSOnly)
         return;
 
+    // Determines the SRS to be used by the geolocation array, if any
+    std::string osGeolocWKT = SRS_WKT_WGS84_LAT_LONG;
+    if (!m_oSRS.IsEmpty())
+    {
+        OGRSpatialReference oGeogCRS;
+        oGeogCRS.CopyGeogCSFrom(&m_oSRS);
+        char *pszWKTTmp = nullptr;
+        const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
+        if (oGeogCRS.exportToWkt(&pszWKTTmp, apszOptions) == OGRERR_NONE)
+        {
+            osGeolocWKT = pszWKTTmp;
+        }
+        CPLFree(pszWKTTmp);
+    }
+
     // Process geolocation arrays from CF "coordinates" attribute.
     std::string osGeolocXName, osGeolocYName;
-    if (ProcessCFGeolocation(nGroupId, nVarId, osGeolocXName, osGeolocYName))
+    if (ProcessCFGeolocation(nGroupId, nVarId, osGeolocWKT, osGeolocXName,
+                             osGeolocYName))
     {
         bool bCanCancelGT = true;
         if ((nVarDimXID != -1) && (nVarDimYID != -1))
@@ -4311,7 +4327,7 @@ void netCDFDataset::SetProjectionFromVar(
         CPLDebug("GDAL_netCDF", "using variables %s and %s for GEOLOCATION",
                  pszGeolocXFullName, pszGeolocYFullName);
 
-        GDALPamDataset::SetMetadataItem("SRS", SRS_WKT_WGS84_LAT_LONG,
+        GDALPamDataset::SetMetadataItem("SRS", osGeolocWKT.c_str(),
                                         "GEOLOCATION");
 
         CPLString osTMP;
@@ -4729,6 +4745,7 @@ bool netCDFDataset::ProcessNASAEMITGeoLocation(int nGroupId, int nVarId)
 }
 
 int netCDFDataset::ProcessCFGeolocation(int nGroupId, int nVarId,
+                                        const std::string &osGeolocWKT,
                                         std::string &osGeolocXNameOut,
                                         std::string &osGeolocYNameOut)
 {
@@ -4801,8 +4818,8 @@ int netCDFDataset::ProcessCFGeolocation(int nGroupId, int nVarId,
                              "using variables %s and %s for GEOLOCATION",
                              pszGeolocXFullName, pszGeolocYFullName);
 
-                    GDALPamDataset::SetMetadataItem(
-                        "SRS", SRS_WKT_WGS84_LAT_LONG, "GEOLOCATION");
+                    GDALPamDataset::SetMetadataItem("SRS", osGeolocWKT.c_str(),
+                                                    "GEOLOCATION");
 
                     CPLString osTMP;
                     osTMP.Printf("NETCDF:\"%s\":%s", osFilename.c_str(),

diff --git a/frmts/netcdf/netcdfdataset.h b/frmts/netcdf/netcdfdataset.h
@@ -502,6 +502,7 @@ class netCDFDataset final : public GDALPamDataset
     bool ProcessNASAEMITGeoLocation(int nGroupId, int nVarId);
 
     int ProcessCFGeolocation(int nGroupId, int nVarId,
+                             const std::string &osGeolocWKT,
                              std::string &osGeolocXNameOut,
                              std::string &osGeolocYNameOut);
     CPLErr Set1DGeolocation(int nGroupId, int nVarId, const char *szDimName);