Merge pull request OSGeo#5098 from rouault/fix_5095

AAIGRID: add support for 'null' as can be generated by D12 software (fixes OSGeo#5095)
rouault · Jan 12, 2022 · e4cdcc7 · e4cdcc7
2 parents 4afab38 + a27636f
commit e4cdcc7
Show file tree

Hide file tree

Showing 2 changed files with 98 additions and 17 deletions.
diff --git a/autotest/gdrivers/aaigrid.py b/autotest/gdrivers/aaigrid.py
@@ -250,7 +250,7 @@ def test_aaigrid_10():
         except OSError:
             pass
 
-    
+
 ###############################################################################
 # Test SIGNIFICANT_DIGITS creation option (same as DECIMAL_PRECISION test)
 
@@ -351,7 +351,52 @@ def test_aaigrid_15():
 
     gdal.Unlink('/vsimem/aaigrid_15.asc')
 
+
+###############################################################################
+# Test support for D12 AAIGRID with null value (#5095)
+
+
+def test_aaigrid_null():
+
+    gdal.FileFromMemBuffer('/vsimem/test_aaigrid_null.asc', """ncols        4
+nrows        1
+xllcorner    0
+yllcorner    -1
+cellsize     1
+NODATA_value  null
+null 1.5 null 3.5
+""")
+
+    ds = gdal.Open('/vsimem/test_aaigrid_null.asc')
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_Float32
+    assert ds.GetRasterBand(1).GetNoDataValue() < -1e38
+    assert ds.GetRasterBand(1).ComputeRasterMinMax() == (1.5,3.5)
+    ds = None
+
+    gdal.Unlink('/vsimem/test_aaigrid_null.asc')
+
+
 ###############################################################################
+# Test support for D12 AAIGRID with null value and force AAIGRID_DATATYPE=Float64 (#5095)
+
 
+def test_aaigrid_null_float64():
+
+    gdal.FileFromMemBuffer('/vsimem/test_aaigrid_null.asc', """ncols        4
+nrows        1
+xllcorner    0
+yllcorner    -1
+cellsize     1
+NODATA_value  null
+null 1.5 null 3.5
+""")
+
+    with gdaltest.config_option('AAIGRID_DATATYPE', 'Float64'):
+        ds = gdal.Open('/vsimem/test_aaigrid_null.asc')
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_Float64
+    assert ds.GetRasterBand(1).GetNoDataValue() < -1e308
+    assert ds.GetRasterBand(1).ComputeRasterMinMax() == (1.5,3.5)
+    ds = None
 
+    gdal.Unlink('/vsimem/test_aaigrid_null.asc')
 
diff --git a/frmts/aaigrid/aaigriddataset.cpp b/frmts/aaigrid/aaigriddataset.cpp
@@ -194,11 +194,25 @@ CPLErr AAIGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
 
         if( pImage != nullptr )
         {
+            // "null" seems to be specific of D12 software
+            // See https://github.com/OSGeo/gdal/issues/5095
             if( eDataType == GDT_Float64 )
-                reinterpret_cast<double *>(pImage)[iPixel] = CPLAtofM(szToken);
+            {
+                if( strcmp(szToken, "null") == 0 )
+                    reinterpret_cast<double *>(pImage)[iPixel] =
+                        -std::numeric_limits<double>::max();
+                else
+                    reinterpret_cast<double *>(pImage)[iPixel] = CPLAtofM(szToken);
+            }
             else if( eDataType == GDT_Float32 )
-                reinterpret_cast<float *>(pImage)[iPixel] =
-                    DoubleToFloatClamp(CPLAtofM(szToken));
+            {
+                if( strcmp(szToken, "null") == 0 )
+                    reinterpret_cast<float *>(pImage)[iPixel] =
+                        -std::numeric_limits<float>::max();
+                else
+                    reinterpret_cast<float *>(pImage)[iPixel] =
+                        DoubleToFloatClamp(CPLAtofM(szToken));
+            }
             else
                 reinterpret_cast<GInt32 *>(pImage)[iPixel] =
                     static_cast<GInt32>(atoi(szToken));
@@ -553,24 +567,42 @@ int AAIGDataset::ParseHeader(const char *pszHeader, const char *pszDataType)
         const char *pszNoData = papszTokens[i + 1];
 
         bNoDataSet = true;
-        dfNoDataValue = CPLAtofM(pszNoData);
-        if( pszDataType == nullptr &&
-            (strchr(pszNoData, '.') != nullptr ||
-             strchr(pszNoData, ',') != nullptr ||
-             std::numeric_limits<int>::min() > dfNoDataValue ||
-             dfNoDataValue > std::numeric_limits<int>::max()) )
+        if( strcmp(pszNoData, "null") == 0 )
         {
-            eDataType = GDT_Float32;
-            if( !CPLIsInf(dfNoDataValue) &&
-                (fabs(dfNoDataValue) < std::numeric_limits<float>::min() ||
-                 fabs(dfNoDataValue) > std::numeric_limits<float>::max()) )
+            // "null" seems to be specific of D12 software
+            // See https://github.com/OSGeo/gdal/issues/5095
+            if( pszDataType == nullptr || eDataType == GDT_Float32 )
             {
+                dfNoDataValue = -std::numeric_limits<float>::max();
+                eDataType = GDT_Float32;
+            }
+            else
+            {
+                dfNoDataValue = -std::numeric_limits<double>::max();
                 eDataType = GDT_Float64;
             }
         }
-        if( eDataType == GDT_Float32 )
+        else
         {
-            dfNoDataValue = MapNoDataToFloat(dfNoDataValue);
+            dfNoDataValue = CPLAtofM(pszNoData);
+            if( pszDataType == nullptr &&
+                (strchr(pszNoData, '.') != nullptr ||
+                 strchr(pszNoData, ',') != nullptr ||
+                 std::numeric_limits<int>::min() > dfNoDataValue ||
+                 dfNoDataValue > std::numeric_limits<int>::max()) )
+            {
+                eDataType = GDT_Float32;
+                if( !CPLIsInf(dfNoDataValue) &&
+                    (fabs(dfNoDataValue) < std::numeric_limits<float>::min() ||
+                     fabs(dfNoDataValue) > std::numeric_limits<float>::max()) )
+                {
+                    eDataType = GDT_Float64;
+                }
+            }
+            if( eDataType == GDT_Float32 )
+            {
+                dfNoDataValue = MapNoDataToFloat(dfNoDataValue);
+            }
         }
     }
 
@@ -1000,7 +1032,11 @@ GDALDataset *AAIGDataset::CommonOpen( GDALOpenInfo *poOpenInfo,
                 poOpenInfo->pabyHeader[i - 1] == '\r' ||
                 poOpenInfo->pabyHeader[i - 2] == '\r' )
             {
-                if( !isalpha(poOpenInfo->pabyHeader[i]) &&
+                if( (!isalpha(poOpenInfo->pabyHeader[i]) ||
+                    // null seems to be specific of D12 software
+                    // See https://github.com/OSGeo/gdal/issues/5095
+                     (i + 5 < poOpenInfo->nHeaderBytes &&
+                      memcmp(poOpenInfo->pabyHeader + i, "null ", 5) == 0)) &&
                     poOpenInfo->pabyHeader[i] != '\n' &&
                     poOpenInfo->pabyHeader[i] != '\r')
                 {