Add GDT_Int64 and GDT_UInt64 data types and handle them in MEM, GTiff…

…, netCDF and Zarr drivers

MIGRATION_GUIDE.TXT

@@ -1,3 +1,9 @@
+MIGRATION GUIDE FROM GDAL 3.4 to GDAL 3.5
+-----------------------------------------
+
+- GDAL drivers may now return raster bands with the new data types GDT_Int64 or
+  GDT_UInt64.
+
 MIGRATION GUIDE FROM GDAL 3.3 to GDAL 3.4
 -----------------------------------------
 

alg/gdal_alg_priv.h

@@ -34,6 +34,8 @@
 
 #ifndef DOXYGEN_SKIP
 
+#include <cstdint>
+
 #include "gdal_alg.h"
 #include "ogr_spatialref.h"
 
@@ -146,10 +148,10 @@ template<class DataType, class EqualityTest> class GDALRasterPolygonEnumeratorT
 
 struct IntEqualityTest
 {
-    bool operator()(GInt32 a, GInt32 b) const { return a == b; }
+    bool operator()(std::int64_t a, std::int64_t b) const { return a == b; }
 };
 
-typedef GDALRasterPolygonEnumeratorT<GInt32, IntEqualityTest> GDALRasterPolygonEnumerator;
+typedef GDALRasterPolygonEnumeratorT<std::int64_t, IntEqualityTest> GDALRasterPolygonEnumerator;
 
 typedef void* (*GDALTransformDeserializeFunc)( CPLXMLNode *psTree );
 

alg/gdalpansharpen.cpp

@@ -267,8 +267,7 @@ GDALPansharpenOperation::Initialize( const GDALPansharpenOptions* psOptionsIn )
     {
         if( psOptionsIn->nBitDepth < 0 || psOptionsIn->nBitDepth > 31 ||
             (eWorkDataType == GDT_Byte && psOptionsIn->nBitDepth > 8) ||
-            (eWorkDataType == GDT_UInt16 && psOptionsIn->nBitDepth > 16) ||
-            (eWorkDataType == GDT_UInt32 && psOptionsIn->nBitDepth > 32) )
+            (eWorkDataType == GDT_UInt16 && psOptionsIn->nBitDepth > 16) )
         {
             CPLError(CE_Failure, CPLE_AppDefined,
                      "Invalid value nBitDepth = %d for type %s",
@@ -282,7 +281,7 @@ GDALPansharpenOperation::Initialize( const GDALPansharpenOptions* psOptionsIn )
         psOptions->nBitDepth = 0;
     if( psOptions->nBitDepth &&
         !(eWorkDataType == GDT_Byte || eWorkDataType == GDT_UInt16 ||
-          eWorkDataType == GDT_UInt32) )
+          eWorkDataType == GDT_UInt32 || eWorkDataType == GDT_UInt64) )
     {
         CPLError(CE_Warning, CPLE_AppDefined,
                  "Ignoring nBitDepth = %d for type %s",
@@ -920,6 +919,18 @@ template<class WorkDataType> CPLErr GDALPansharpenOperation::WeightedBrovey(
                            nValues, nBandValues, nMaxValue);
             break;
 
+        case GDT_UInt64:
+            WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
+                           static_cast<std::uint64_t *>(pDataBuf),
+                           nValues, nBandValues, nMaxValue);
+            break;
+
+        case GDT_Int64:
+            WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
+                           static_cast<std::int64_t *>(pDataBuf),
+                           nValues, nBandValues, nMaxValue);
+            break;
+
         case GDT_Float32:
             WeightedBrovey(pPanBuffer, pUpsampledSpectralBuffer,
                            static_cast<float *>(pDataBuf),
@@ -983,6 +994,18 @@ template<class WorkDataType> CPLErr GDALPansharpenOperation::WeightedBrovey(
                 static_cast<GInt32 *>(pDataBuf), nValues, nBandValues, 0);
             break;
 
+        case GDT_UInt64:
+            WeightedBrovey3<WorkDataType, std::uint64_t, FALSE>(
+                pPanBuffer, pUpsampledSpectralBuffer,
+                static_cast<std::uint64_t *>(pDataBuf), nValues, nBandValues, 0);
+            break;
+
+        case GDT_Int64:
+            WeightedBrovey3<WorkDataType, std::int64_t, FALSE>(
+                pPanBuffer, pUpsampledSpectralBuffer,
+                static_cast<std::int64_t *>(pDataBuf), nValues, nBandValues, 0);
+            break;
+
         case GDT_Float32:
             WeightedBrovey3<WorkDataType, float, FALSE>(
                 pPanBuffer, pUpsampledSpectralBuffer,
@@ -1699,6 +1722,20 @@ GDALPansharpenOperation::PansharpenChunk( GDALDataType eWorkDataType,
                                   nValues, nBandValues);
             break;
 
+        case GDT_UInt64:
+            eErr = WeightedBrovey(static_cast<const std::uint64_t*>(pPanBuffer),
+                                  static_cast<const std::uint64_t*>(pUpsampledSpectralBuffer),
+                                  pDataBuf, eBufDataType,
+                                  nValues, nBandValues, nMaxValue);
+            break;
+
+        case GDT_Int64:
+            eErr = WeightedBrovey(static_cast<const std::int64_t*>(pPanBuffer),
+                                  static_cast<const std::int64_t*>(pUpsampledSpectralBuffer),
+                                  pDataBuf, eBufDataType,
+                                  nValues, nBandValues);
+            break;
+
         case GDT_Float32:
             eErr = WeightedBrovey(static_cast<const float*>(pPanBuffer),
                                   static_cast<const float*>(pUpsampledSpectralBuffer),

alg/gdalrasterize.cpp

@@ -137,6 +137,12 @@ void gvBurnScanline( void *pCBData, int nY, int nXStart, int nXEnd,
         case GDT_UInt32:
             gvBurnScanlineBasic<GUInt32>( psInfo, nY, nXStart, nXEnd, dfVariant );
             break;
+        case GDT_Int64:
+            gvBurnScanlineBasic<std::int64_t>( psInfo, nY, nXStart, nXEnd, dfVariant );
+            break;
+        case GDT_UInt64:
+            gvBurnScanlineBasic<std::uint64_t>( psInfo, nY, nXStart, nXEnd, dfVariant );
+            break;
         case GDT_Float32:
             gvBurnScanlineBasic<float>( psInfo, nY, nXStart, nXEnd, dfVariant );
             break;
@@ -168,7 +174,8 @@ void gvBurnPointBasic( GDALRasterizeInfo *psInfo,
                                  + nY * psInfo->nLineSpace + nX * psInfo->nPixelSpace;
 
         T* pbyPixel = reinterpret_cast<T*>(pbyInsert);
-        burnValue += ( psInfo->eMergeAlg != GRMA_Add ) ? 0 : *pbyPixel;
+        if( psInfo->eMergeAlg == GRMA_Add )
+            burnValue += static_cast<double>(*pbyPixel);
         GDALCopyWord(burnValue, *pbyPixel);
     }
 }
@@ -203,6 +210,12 @@ void gvBurnPoint( void *pCBData, int nY, int nX, double dfVariant )
         case GDT_UInt32:
             gvBurnPointBasic<GUInt32>( psInfo, nY, nX, dfVariant );
             break;
+        case GDT_Int64:
+            gvBurnPointBasic<std::int64_t>( psInfo, nY, nX, dfVariant );
+            break;
+        case GDT_UInt64:
+            gvBurnPointBasic<std::uint64_t>( psInfo, nY, nX, dfVariant );
+            break;
         case GDT_Float32:
             gvBurnPointBasic<float>( psInfo, nY, nX, dfVariant );
             break;

alg/gdalrasterpolygonenumerator.cpp

@@ -299,7 +299,7 @@ void GDALRasterPolygonEnumeratorT<DataType, EqualityTest>::ProcessLine(
     }
 }
 
-template class GDALRasterPolygonEnumeratorT<GInt32, IntEqualityTest>;
+template class GDALRasterPolygonEnumeratorT<std::int64_t, IntEqualityTest>;
 
 template class GDALRasterPolygonEnumeratorT<float, FloatEqualityTest>;
 

alg/gdalsievefilter.cpp

@@ -77,7 +77,7 @@ CPL_CVSID("$Id$")
 static CPLErr
 GPMaskImageData( GDALRasterBandH hMaskBand, GByte *pabyMaskLine,
                  int iY, int nXSize,
-                 GInt32 *panImageLine )
+                 std::int64_t *panImageLine )
 
 {
     const CPLErr eErr =
@@ -113,7 +113,7 @@ GPMaskImageData( GDALRasterBandH hMaskBand, GByte *pabyMaskLine,
 
 static inline void CompareNeighbour( int nPolyId1, int nPolyId2,
                                      int *panPolyIdMap,
-                                     int * /* panPolyValue */,
+                                     std::int64_t * /* panPolyValue */,
                                      std::vector<int> &anPolySizes,
                                      std::vector<int> &anBigNeighbour )
 
@@ -215,16 +215,16 @@ GDALSieveFilter( GDALRasterBandH hSrcBand, GDALRasterBandH hMaskBand,
 /* -------------------------------------------------------------------- */
     int nXSize = GDALGetRasterBandXSize( hSrcBand );
     int nYSize = GDALGetRasterBandYSize( hSrcBand );
-    GInt32 *panLastLineVal = static_cast<GInt32 *>(
+    auto *panLastLineVal = static_cast<std::int64_t *>(
+        VSI_MALLOC2_VERBOSE(sizeof(std::int64_t), nXSize));
+    auto *panThisLineVal = static_cast<std::int64_t *>(
+        VSI_MALLOC2_VERBOSE(sizeof(std::int64_t), nXSize));
+    auto *panLastLineId = static_cast<GInt32 *>(
         VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
-    GInt32 *panThisLineVal = static_cast<GInt32 *>(
-        VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
-    GInt32 *panLastLineId = static_cast<GInt32 *>(
-        VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
-    GInt32 *panThisLineId = static_cast<GInt32 *>(
-        VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
-    GInt32 *panThisLineWriteVal = static_cast<GInt32 *>(
+    auto *panThisLineId = static_cast<GInt32 *>(
         VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize));
+    auto *panThisLineWriteVal = static_cast<std::int64_t *>(
+        VSI_MALLOC2_VERBOSE(sizeof(std::int64_t), nXSize));
     GByte *pabyMaskLine =
         hMaskBand != nullptr
         ? static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize))
@@ -257,7 +257,7 @@ GDALSieveFilter( GDALRasterBandH hSrcBand, GDALRasterBandH hMaskBand,
         eErr = GDALRasterIO(
             hSrcBand,
             GF_Read, 0, iY, nXSize, 1,
-            panThisLineVal, nXSize, 1, GDT_Int32, 0, 0 );
+            panThisLineVal, nXSize, 1, GDT_Int64, 0, 0 );
 
         if( eErr == CE_None && hMaskBand != nullptr )
             eErr = GPMaskImageData(hMaskBand, pabyMaskLine, iY, nXSize,
@@ -355,7 +355,7 @@ GDALSieveFilter( GDALRasterBandH hSrcBand, GDALRasterBandH hMaskBand,
 /*      Read the image data.                                            */
 /* -------------------------------------------------------------------- */
         eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1,
-                             panThisLineVal, nXSize, 1, GDT_Int32, 0, 0 );
+                             panThisLineVal, nXSize, 1, GDT_Int64, 0, 0 );
 
         if( eErr == CE_None && hMaskBand != nullptr )
             eErr = GPMaskImageData( hMaskBand, pabyMaskLine, iY, nXSize,
@@ -532,9 +532,9 @@ GDALSieveFilter( GDALRasterBandH hSrcBand, GDALRasterBandH hMaskBand,
 /*      Read the image data.                                            */
 /* -------------------------------------------------------------------- */
         eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1,
-                             panThisLineVal, nXSize, 1, GDT_Int32, 0, 0 );
+                             panThisLineVal, nXSize, 1, GDT_Int64, 0, 0 );
 
-        memcpy( panThisLineWriteVal, panThisLineVal, 4 * nXSize );
+        memcpy( panThisLineWriteVal, panThisLineVal, sizeof(panThisLineVal[0]) * nXSize );
 
         if( eErr == CE_None && hMaskBand != nullptr )
             eErr = GPMaskImageData( hMaskBand, pabyMaskLine, iY, nXSize,
@@ -580,7 +580,7 @@ GDALSieveFilter( GDALRasterBandH hSrcBand, GDALRasterBandH hMaskBand,
 /*      Write the update data out.                                      */
 /* -------------------------------------------------------------------- */
         eErr = GDALRasterIO( hDstBand, GF_Write, 0, iY, nXSize, 1,
-                             panThisLineWriteVal, nXSize, 1, GDT_Int32, 0, 0 );
+                             panThisLineWriteVal, nXSize, 1, GDT_Int64, 0, 0 );
 
 /* -------------------------------------------------------------------- */
 /*      Swap pixel value, and polygon id lines to be ready for the      */

alg/gdalwarpkernel.cpp

@@ -1591,6 +1591,16 @@ static bool GWKSetPixelValue( const GDALWarpKernel *poWK, int iBand,
             dfDstImag = 0.0;
             break;
 
+          case GDT_Int64:
+            dfDstReal = static_cast<double>(reinterpret_cast<std::int64_t*>(pabyDst)[iDstOffset]);
+            dfDstImag = 0.0;
+            break;
+
+          case GDT_UInt64:
+            dfDstReal = static_cast<double>(reinterpret_cast<std::uint64_t*>(pabyDst)[iDstOffset]);
+            dfDstImag = 0.0;
+            break;
+
           case GDT_Float32:
             dfDstReal = reinterpret_cast<float*>(pabyDst)[iDstOffset];
             dfDstImag = 0.0;
@@ -1690,6 +1700,14 @@ static bool GWKSetPixelValue( const GDALWarpKernel *poWK, int iBand,
         CLAMP(GInt32);
         break;
 
+      case GDT_UInt64:
+        CLAMP(std::uint64_t);
+        break;
+
+      case GDT_Int64:
+        CLAMP(std::int64_t);
+        break;
+
       case GDT_Float32:
         reinterpret_cast<float*>(pabyDst)[iDstOffset] = static_cast<float>(dfReal);
         break;
@@ -1810,6 +1828,14 @@ static bool GWKSetPixelValueReal( const GDALWarpKernel *poWK, int iBand,
             dfDstReal = reinterpret_cast<GUInt32*>(pabyDst)[iDstOffset];
             break;
 
+          case GDT_Int64:
+            dfDstReal = static_cast<double>(reinterpret_cast<std::int64_t*>(pabyDst)[iDstOffset]);
+            break;
+
+          case GDT_UInt64:
+            dfDstReal = static_cast<double>(reinterpret_cast<std::uint64_t*>(pabyDst)[iDstOffset]);
+            break;
+
           case GDT_Float32:
             dfDstReal = reinterpret_cast<float*>(pabyDst)[iDstOffset];
             break;
@@ -1861,6 +1887,14 @@ static bool GWKSetPixelValueReal( const GDALWarpKernel *poWK, int iBand,
         CLAMP(GInt32);
         break;
 
+      case GDT_UInt64:
+        CLAMP(std::uint64_t);
+        break;
+
+      case GDT_Int64:
+        CLAMP(std::int64_t);
+        break;
+
       case GDT_Float32:
         reinterpret_cast<float*>(pabyDst)[iDstOffset] = static_cast<float>(dfReal);
         break;
@@ -1926,6 +1960,16 @@ static bool GWKGetPixelValue( const GDALWarpKernel *poWK, int iBand,
         *pdfImag = 0.0;
         break;
 
+      case GDT_Int64:
+        *pdfReal = static_cast<double>(reinterpret_cast<std::int64_t*>(pabySrc)[iSrcOffset]);
+        *pdfImag = 0.0;
+        break;
+
+      case GDT_UInt64:
+        *pdfReal = static_cast<double>(reinterpret_cast<std::uint64_t*>(pabySrc)[iSrcOffset]);
+        *pdfImag = 0.0;
+        break;
+
       case GDT_Float32:
         *pdfReal = reinterpret_cast<float*>(pabySrc)[iSrcOffset];
         *pdfImag = 0.0;
@@ -2011,6 +2055,14 @@ static bool GWKGetPixelValueReal( const GDALWarpKernel *poWK, int iBand,
         *pdfReal = reinterpret_cast<GUInt32*>(pabySrc)[iSrcOffset];
         break;
 
+      case GDT_Int64:
+        *pdfReal = static_cast<double>(reinterpret_cast<std::int64_t*>(pabySrc)[iSrcOffset]);
+        break;
+
+      case GDT_UInt64:
+        *pdfReal = static_cast<double>(reinterpret_cast<std::uint64_t*>(pabySrc)[iSrcOffset]);
+        break;
+
       case GDT_Float32:
         *pdfReal = reinterpret_cast<float*>(pabySrc)[iSrcOffset];
         break;
@@ -2172,6 +2224,30 @@ static bool GWKGetPixelRow( const GDALWarpKernel *poWK, int iBand,
             break;
         }
 
+        case GDT_Int64:
+        {
+            auto pSrc = reinterpret_cast<std::int64_t*>(poWK->papabySrcImage[iBand]);
+            pSrc += iSrcOffset;
+            for( int i = 0; i < nSrcLen; i += 2 )
+            {
+                adfReal[i] = static_cast<double>(pSrc[i]);
+                adfReal[i+1] = static_cast<double>(pSrc[i+1]);
+            }
+            break;
+        }
+
+        case GDT_UInt64:
+        {
+            auto pSrc = reinterpret_cast<std::uint64_t*>(poWK->papabySrcImage[iBand]);
+            pSrc += iSrcOffset;
+            for( int i = 0; i < nSrcLen; i += 2 )
+            {
+                adfReal[i] = static_cast<double>(pSrc[i]);
+                adfReal[i+1] = static_cast<double>(pSrc[i+1]);
+            }
+            break;
+        }
+
         case GDT_Float32:
         {
             float* pSrc = reinterpret_cast<float*>(poWK->papabySrcImage[iBand]);