netCDF: fix reading MODIS_ARRAY.nc file that has non CF-compliant way…

… of expressing axis and CRS (fixes #4075)
OSGeo · Jul 7, 2021 · 080e975 · 080e975
1 parent a26ec1e
commit 080e975
Show file tree

Hide file tree

Showing 3 changed files with 39 additions and 2 deletions.
diff --git a/autotest/gdrivers/data/netcdf/MODIS_ARRAY.nc b/autotest/gdrivers/data/netcdf/MODIS_ARRAY.nc
diff --git a/autotest/gdrivers/netcdf.py b/autotest/gdrivers/netcdf.py
@@ -4687,6 +4687,17 @@ def test_netcdf_metadata_sentinel5():
     assert j['attributes']['ISO_METADATA'] == expected
 
 
+###############################################################################
+# Test opening a file with particular georeferencing encoding
+
+
+def test_netcdf_modis_array():
+
+    ds = gdal.Open('data/netcdf/MODIS_ARRAY.nc')
+    assert ds.GetGeoTransform(can_return_null=True) is not None
+    assert ds.GetSpatialRef() is not None
+
+
 def test_clean_tmp():
     # [KEEP THIS AS THE LAST TEST]
     # i.e. please do not add any tests after this one. Put new ones above.

diff --git a/gdal/frmts/netcdf/netcdfdataset.cpp b/gdal/frmts/netcdf/netcdfdataset.cpp
@@ -3368,6 +3368,17 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
             // oSRS.SetWellKnownGeogCS("WGS84");
         }
     }
+    else
+    {
+        // Dataset from https://github.com/OSGeo/gdal/issues/4075 has a "crs"
+        // attribute hold on the variable of interest that contains a PROJ.4 string
+        pszValue = FetchAttr(nGroupId, nVarId, "crs");
+        if( pszValue && strstr(pszValue, "+proj=") &&
+            oSRS.importFromProj4(pszValue) == OGRERR_NONE )
+        {
+            bGotCfSRS = true;
+        }
+    }
     // Read projection coordinates.
 
     int nGroupDimXID = -1;
@@ -3392,6 +3403,11 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
         // variables without same name than dimension using the same resolving
         // logic. This should handle for example NASA Ocean Color L2 products.
 
+        const bool bIgnoreXYAxisNameChecks =
+            CPLTestBool(CPLGetConfigOption("GDAL_NETCDF_IGNORE_XY_AXIS_NAME_CHECKS", "NO")) ||
+            // Dataset from https://github.com/OSGeo/gdal/issues/4075 has a res and transform attributes
+            (FetchAttr(nGroupId, nVarId, "res") != nullptr &&
+             FetchAttr(nGroupId, nVarId, "transform") != nullptr);
 
         // Check that they are 1D or 2D variables
         if( nVarDimXID >= 0 )
@@ -3400,14 +3416,19 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
             nc_inq_varndims(nGroupId, nVarDimXID, &ndims);
             if( ndims == 0 || ndims > 2 )
                 nVarDimXID = -1;
-            else
+            else if( !bIgnoreXYAxisNameChecks )
             {
                 if( !NCDFIsVarLongitude(nGroupId, nVarDimXID, nullptr) &&
                     !NCDFIsVarProjectionX(nGroupId, nVarDimXID, nullptr) &&
                     // In case of inversion of X/Y
                     !NCDFIsVarLatitude(nGroupId, nVarDimXID, nullptr) &&
                     !NCDFIsVarProjectionY(nGroupId, nVarDimXID, nullptr) )
                 {
+                    CPLDebug("netCDF",
+                             "Georeferencing ignored due to non-specific "
+                             "enough X axis name. "
+                             "Set GDAL_NETCDF_IGNORE_XY_AXIS_NAME_CHECKS=YES "
+                             "as configuration option to bypass this check");
                     nVarDimXID = -1;
                 }
             }
@@ -3419,14 +3440,19 @@ void netCDFDataset::SetProjectionFromVar( int nGroupId, int nVarId,
             nc_inq_varndims(nGroupId, nVarDimYID, &ndims);
             if( ndims == 0 || ndims > 2 )
                 nVarDimYID = -1;
-            else
+            else if( !bIgnoreXYAxisNameChecks )
             {
                 if( !NCDFIsVarLatitude(nGroupId, nVarDimYID, nullptr) &&
                     !NCDFIsVarProjectionY(nGroupId, nVarDimYID, nullptr) &&
                     // In case of inversion of X/Y
                     !NCDFIsVarLongitude(nGroupId, nVarDimYID, nullptr) &&
                     !NCDFIsVarProjectionX(nGroupId, nVarDimYID, nullptr) )
                 {
+                    CPLDebug("netCDF",
+                             "Georeferencing ignored due to non-specific "
+                             "enough Y axis name. "
+                             "Set GDAL_NETCDF_IGNORE_XY_AXIS_NAME_CHECKS=YES "
+                             "as configuration option to bypass this check");
                     nVarDimYID = -1;
                 }
             }