update to GDAL 3.8.5 (#147)

* update GDAL to 3.8.5

* add flatbbuffers define

* add missing source file

CHANGELOG.md

@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [3.8.5] 2024-04-09
+
+### Added
+ - GDAL 3.8.5
+
+
 ## [3.8.4] 2024-02-19
 
 ### Added

deps/libgdal.sh

@@ -5,7 +5,7 @@ set -eu
 DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
 cd "$DIR/libgdal"
 
-GDAL_VERSION=3.8.4
+GDAL_VERSION=3.8.5
 GDAL_VERSION_SUFFIX=
 dir_gdal=./gdal
 dir_formats_gyp=./gyp-formats

deps/libgdal/common.gypi

@@ -47,7 +47,8 @@
 			"ZIP_SUPPORT=1",
       "HAVE_LIBZ=1",
 			"JPEG_SUPPORT=1",
-      "RENAME_INTERNAL_SHAPELIB_SYMBOLS=1"
+      "RENAME_INTERNAL_SHAPELIB_SYMBOLS=1",
+      "flatbuffers=gdal_flatbuffers"
 		],
 		"dependencies": [
 			"<(deps_dir)/libexpat/libexpat.gyp:libexpat",

deps/libgdal/gdal/CITATION.cff

@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: Please cite this software using these metadata or in the CITATION file.
 type: software
 title: GDAL
-version: 3.8.4
-date-released: 2024-02-08
+version: 3.8.5
+date-released: 2024-04-02
 doi: 10.5281/zenodo.5884351
 abstract: GDAL is a translator library for raster and vector geospatial data
   formats that is released under an MIT style Open Source License by the Open

deps/libgdal/gdal/CMakeLists.txt

@@ -174,6 +174,16 @@ set_package_properties(
   URL "http://swig.org/"
   TYPE RECOMMENDED)
 
+# Workaround Intel MKL library outputting a message in stdout, which cause
+# incorrect detection of numpy.get_include() and numpy.__version__
+# See https://github.com/numpy/numpy/issues/23775
+# Fix also proposed to upstream CMake FindPython per
+# https://gitlab.kitware.com/cmake/cmake/-/merge_requests/9359/
+if(DEFINED ENV{MKL_ENABLE_INSTRUCTIONS})
+  set(BACKUP_ENV_VAR_MKL_ENABLE_INSTRUCTIONS ENV{MKL_ENABLE_INSTRUCTIONS})
+endif()
+set(ENV{MKL_ENABLE_INSTRUCTIONS} "SSE4_2")
+
 if (Python_LOOKUP_VERSION)
     set(Python_FIND_STRATEGY VERSION)
     find_package(Python ${Python_LOOKUP_VERSION} EXACT COMPONENTS Interpreter Development NumPy)
@@ -227,6 +237,14 @@ if (Python_Interpreter_FOUND)
   endif()
 endif()
 
+# Restore previous value of MKL_ENABLE_INSTRUCTIONS
+if(DEFINED BACKUP_ENV_VAR_MKL_ENABLE_INSTRUCTIONS)
+  set(ENV{MKL_ENABLE_INSTRUCTIONS} ${BACKUP_ENV_VAR_MKL_ENABLE_INSTRUCTIONS})
+else()
+  unset(ENV{MKL_ENABLE_INSTRUCTIONS})
+endif()
+
+
 if (SWIG_FOUND AND (Python_Interpreter_FOUND OR Python_FOUND))
   option(BUILD_PYTHON_BINDINGS "Build Python bindings" ON)
 else()

deps/libgdal/gdal/VERSION

@@ -1 +1 @@
-3.8.4
+3.8.5

deps/libgdal/gdal/alg/gdalpansharpen.cpp

@@ -308,7 +308,7 @@ GDALPansharpenOperation::Initialize(const GDALPansharpenOptions *psOptionsIn)
             CPLError(CE_Failure, CPLE_AppDefined,
                      "Dimensions of input spectral band %d different from "
                      "first spectral band",
-                     i);
+                     i + 1);
             return CE_Failure;
         }
 
@@ -318,8 +318,8 @@ GDALPansharpenOperation::Initialize(const GDALPansharpenOptions *psOptionsIn)
         {
             CPLError(CE_Failure, CPLE_AppDefined,
                      "Cannot retrieve dataset associated with "
-                     "ahInputSpectralBands[%i]",
-                     i);
+                     "input spectral band %d",
+                     i + 1);
             return CE_Failure;
         }
         // Make sure that the band is really a first level child of the owning dataset
@@ -335,16 +335,16 @@ GDALPansharpenOperation::Initialize(const GDALPansharpenOptions *psOptionsIn)
         {
             CPLError(
                 CE_Failure, CPLE_AppDefined,
-                "ahInputSpectralBands[%d] band has no associated geotransform",
-                i);
+                "input spectral band %d band has no associated geotransform",
+                i + 1);
             return CE_Failure;
         }
         if (adfMSGT != adfRefMSGT)
         {
             CPLError(CE_Failure, CPLE_AppDefined,
-                     "ahInputSpectralBands[%d] band has a different "
-                     "geotransform than ahInputSpectralBands[0]",
-                     i);
+                     "input spectral band %d has a different "
+                     "geotransform than the first spectral band",
+                     i + 1);
             return CE_Failure;
         }
 

deps/libgdal/gdal/apps/gdalinfo_lib.cpp

@@ -341,19 +341,33 @@ char *GDALInfo(GDALDatasetH hDataset, const GDALInfoOptions *psOptions)
 
     if (bJson)
     {
-        json_object *poSize = json_object_new_array();
-        json_object *poSizeX =
-            json_object_new_int(GDALGetRasterXSize(hDataset));
-        json_object *poSizeY =
-            json_object_new_int(GDALGetRasterYSize(hDataset));
-
-        json_object_array_add(poSize, poSizeX);
-        json_object_array_add(poSize, poSizeY);
-
-        json_object *poStacSize = nullptr;
-        json_object_deep_copy(poSize, &poStacSize, nullptr);
-        json_object_object_add(poJsonObject, "size", poSize);
-        json_object_object_add(poStac, "proj:shape", poStacSize);
+        {
+            json_object *poSize = json_object_new_array();
+            json_object *poSizeX =
+                json_object_new_int(GDALGetRasterXSize(hDataset));
+            json_object *poSizeY =
+                json_object_new_int(GDALGetRasterYSize(hDataset));
+
+            // size is X, Y ordered
+            json_object_array_add(poSize, poSizeX);
+            json_object_array_add(poSize, poSizeY);
+
+            json_object_object_add(poJsonObject, "size", poSize);
+        }
+
+        {
+            json_object *poStacSize = json_object_new_array();
+            json_object *poSizeX =
+                json_object_new_int(GDALGetRasterXSize(hDataset));
+            json_object *poSizeY =
+                json_object_new_int(GDALGetRasterYSize(hDataset));
+
+            // ... but ... proj:shape is Y, X ordered.
+            json_object_array_add(poStacSize, poSizeY);
+            json_object_array_add(poStacSize, poSizeX);
+
+            json_object_object_add(poStac, "proj:shape", poStacSize);
+        }
     }
     else
     {
@@ -702,15 +716,20 @@ char *GDALInfo(GDALDatasetH hDataset, const GDALInfoOptions *psOptions)
     {
         OGRSpatialReferenceH hLatLong = nullptr;
 
-        OGRErr eErr = OGRERR_NONE;
-        // Check that it looks like Earth before trying to reproject to wgs84...
-        if (bJson &&
-            fabs(OSRGetSemiMajor(hProj, &eErr) - 6378137.0) < 10000.0 &&
-            eErr == OGRERR_NONE)
+        if (bJson)
         {
-            bTransformToWGS84 = true;
-            hLatLong = OSRNewSpatialReference(nullptr);
-            OSRSetWellKnownGeogCS(hLatLong, "WGS84");
+            // Check that it looks like Earth before trying to reproject to wgs84...
+            // OSRGetSemiMajor() may raise an error on CRS like Engineering CRS
+            CPLErrorHandlerPusher oPusher(CPLQuietErrorHandler);
+            CPLErrorStateBackuper oCPLErrorHandlerPusher;
+            OGRErr eErr = OGRERR_NONE;
+            if (fabs(OSRGetSemiMajor(hProj, &eErr) - 6378137.0) < 10000.0 &&
+                eErr == OGRERR_NONE)
+            {
+                bTransformToWGS84 = true;
+                hLatLong = OSRNewSpatialReference(nullptr);
+                OSRSetWellKnownGeogCS(hLatLong, "WGS84");
+            }
         }
         else
         {

deps/libgdal/gdal/apps/gdalwarp_lib.cpp

@@ -4790,16 +4790,56 @@ static void RemoveZeroWidthSlivers(OGRGeometry *poGeom)
     const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
     if (eType == wkbMultiPolygon)
     {
-        for (auto poSubGeom : *(poGeom->toMultiPolygon()))
+        auto poMP = poGeom->toMultiPolygon();
+        int nNumGeometries = poMP->getNumGeometries();
+        for (int i = 0; i < nNumGeometries; /* incremented in loop */)
         {
-            RemoveZeroWidthSlivers(poSubGeom);
+            auto poPoly = poMP->getGeometryRef(i);
+            RemoveZeroWidthSlivers(poPoly);
+            if (poPoly->IsEmpty())
+            {
+                CPLDebug("WARP",
+                         "RemoveZeroWidthSlivers: removing empty polygon");
+                poMP->removeGeometry(i, /* bDelete = */ true);
+                --nNumGeometries;
+            }
+            else
+            {
+                ++i;
+            }
         }
     }
     else if (eType == wkbPolygon)
     {
-        for (auto poSubGeom : *(poGeom->toPolygon()))
+        auto poPoly = poGeom->toPolygon();
+        if (auto poExteriorRing = poPoly->getExteriorRing())
         {
-            RemoveZeroWidthSlivers(poSubGeom);
+            RemoveZeroWidthSlivers(poExteriorRing);
+            if (poExteriorRing->getNumPoints() < 4)
+            {
+                poPoly->empty();
+                return;
+            }
+        }
+        int nNumInteriorRings = poPoly->getNumInteriorRings();
+        for (int i = 0; i < nNumInteriorRings; /* incremented in loop */)
+        {
+            auto poRing = poPoly->getInteriorRing(i);
+            RemoveZeroWidthSlivers(poRing);
+            if (poRing->getNumPoints() < 4)
+            {
+                CPLDebug(
+                    "WARP",
+                    "RemoveZeroWidthSlivers: removing empty interior ring");
+                constexpr int OFFSET_EXTERIOR_RING = 1;
+                poPoly->removeRing(i + OFFSET_EXTERIOR_RING,
+                                   /* bDelete = */ true);
+                --nNumInteriorRings;
+            }
+            else
+            {
+                ++i;
+            }
         }
     }
     else if (eType == wkbLineString)

deps/libgdal/gdal/cmake/helpers/configure.cmake

@@ -257,7 +257,7 @@ else ()
   # available at build time.
   # This is also necessary for Mac Catalyst builds.
   # Cf https://lists.osgeo.org/pipermail/gdal-dev/2022-August/056174.html
-  if (${CMAKE_SYSTEM_NAME} MATCHES "iOS" OR ${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+  if (${CMAKE_SYSTEM_NAME} MATCHES "iOS|visionOS|tvOS|watchOS|Darwin")
     set(VSI_FOPEN64 "fopen")
     set(VSI_FTRUNCATE64 "ftruncate")
     set(VSI_FTELL64 "ftell")

deps/libgdal/gdal/cmake/modules/packages/FindGEOS.cmake

@@ -19,8 +19,8 @@
 
 find_program(GEOS_CONFIG geos-config)
 if(GEOS_CONFIG)
-    exec_program(${GEOS_CONFIG} ARGS --version OUTPUT_VARIABLE GEOS_VERSION)
-    exec_program(${GEOS_CONFIG} ARGS --prefix OUTPUT_VARIABLE GEOS_PREFIX)
+    execute_process(COMMAND "${GEOS_CONFIG}" --version OUTPUT_VARIABLE GEOS_VERSION)
+    execute_process(COMMAND "${GEOS_CONFIG}" --prefix OUTPUT_VARIABLE GEOS_PREFIX)
 endif()
 
 find_path(GEOS_INCLUDE_DIR NAMES geos_c.h

deps/libgdal/gdal/data/gdalinfo_output.schema.json

@@ -176,6 +176,7 @@
           }
         },
         "size": {
+          "$comment": "note that the order of items in side is width,height",
           "$ref": "#/definitions/arrayOfTwoIntegers"
         },
         "coordinateSystem": {
@@ -306,6 +307,7 @@
         },
 
         "proj:shape": {
+          "$comment": "note that the order of items in proj:shape is height,width starting with GDAL 3.8.5 (previous versions ordered it wrongly as width,height)",
           "title": "Shape",
           "type": "array",
           "minItems": 2,

deps/libgdal/gdal/data/gdalvrt.xsd

@@ -138,13 +138,15 @@
     <xs:complexType name="PanchroBandType">
         <xs:sequence>
             <xs:element name="SourceFilename" type="SourceFilenameType"/>
+            <xs:element name="OpenOptions" type="OpenOptionsType"/>
             <xs:element name="SourceBand" type="xs:string"/>  <!-- should be refined into xs:nonNegativeInteger or mask,xs:nonNegativeInteger -->
         </xs:sequence>
     </xs:complexType>
 
     <xs:complexType name="SpectralBandType">
         <xs:sequence>
             <xs:element name="SourceFilename" type="SourceFilenameType"/>
+            <xs:element name="OpenOptions" type="OpenOptionsType"/>
             <xs:element name="SourceBand" type="xs:string"/>  <!-- should be refined into xs:nonNegativeInteger or mask,xs:nonNegativeInteger -->
         </xs:sequence>
         <xs:attribute name="dstBand" type="xs:nonNegativeInteger"/>

deps/libgdal/gdal/frmts/dimap/dimapdataset.cpp

@@ -1515,6 +1515,8 @@ int DIMAPDataset::ReadImageInformation2()
         "",
         "Processing_Information.Product_Settings.Geometric_Settings",
         "GEOMETRIC_",
+        "Processing_Information.Product_Settings.Radiometric_Settings",
+        "RADIOMETRIC_",
         "Quality_Assessment.Imaging_Quality_Measurement",
         "CLOUDCOVER_",
         nullptr,

deps/libgdal/gdal/frmts/ers/ersdataset.cpp

@@ -881,17 +881,18 @@ int ERSProxyRasterBand::GetOverviewCount()
 GDALDataset *ERSDataset::Open(GDALOpenInfo *poOpenInfo)
 
 {
+    if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
+        return nullptr;
+
+    int &nRecLevel = GetRecLevel();
     // cppcheck-suppress knownConditionTrueFalse
-    if (GetRecLevel())
+    if (nRecLevel)
     {
         CPLError(CE_Failure, CPLE_AppDefined,
                  "Attempt at recursively opening ERS dataset");
         return nullptr;
     }
 
-    if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
-        return nullptr;
-
     /* -------------------------------------------------------------------- */
     /*      Ingest the file as a tree of header nodes.                      */
     /* -------------------------------------------------------------------- */
@@ -1020,7 +1021,6 @@ GDALDataset *ERSDataset::Open(GDALOpenInfo *poOpenInfo)
     /* -------------------------------------------------------------------- */
     if (EQUAL(poHeader->Find("DataSetType", ""), "Translated"))
     {
-        int &nRecLevel = GetRecLevel();
         nRecLevel++;
         poDS->poDepFile = GDALDataset::FromHandle(
             GDALOpen(osDataFilePath, poOpenInfo->eAccess));

deps/libgdal/gdal/frmts/grib/gribdataset.cpp

@@ -1760,13 +1760,15 @@ void GRIBArray::Init(GRIBGroup *poGroup, GRIBDataset *poDS,
         {
             bool bOK = true;
             auto poVar = oIterX->second->GetIndexingVariable();
+            constexpr double EPSILON = 1e-10;
             if (poVar)
             {
                 GUInt64 nStart = 0;
                 size_t nCount = 1;
                 double dfVal = 0;
                 poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt, &dfVal);
-                if (dfVal != adfGT[0] + 0.5 * adfGT[1])
+                if (std::fabs(dfVal - (adfGT[0] + 0.5 * adfGT[1])) >
+                    EPSILON * std::fabs(dfVal))
                 {
                     bOK = false;
                 }
@@ -1781,8 +1783,10 @@ void GRIBArray::Init(GRIBGroup *poGroup, GRIBDataset *poDS,
                     double dfVal = 0;
                     poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt,
                                 &dfVal);
-                    if (dfVal != adfGT[3] + poDS->nRasterYSize * adfGT[5] -
-                                     0.5 * adfGT[5])
+                    if (std::fabs(dfVal -
+                                  (adfGT[3] + poDS->nRasterYSize * adfGT[5] -
+                                   0.5 * adfGT[5])) >
+                        EPSILON * std::fabs(dfVal))
                     {
                         bOK = false;
                     }