diff --git a/MIGRATION_GUIDE.TXT b/MIGRATION_GUIDE.TXT
index a37f4d8c2c17..85870aba38d4 100644
--- a/MIGRATION_GUIDE.TXT
+++ b/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
 -----------------------------------------
 
diff --git a/alg/gdal_alg_priv.h b/alg/gdal_alg_priv.h
index 3d97c3345ba2..1eee84f1ecbf 100644
--- a/alg/gdal_alg_priv.h
+++ b/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 );
 
diff --git a/alg/gdalpansharpen.cpp b/alg/gdalpansharpen.cpp
index cb30f895c87b..6029142e8556 100644
--- a/alg/gdalpansharpen.cpp
+++ b/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),
diff --git a/alg/gdalrasterize.cpp b/alg/gdalrasterize.cpp
index a5ee85033d94..8d3935475603 100644
--- a/alg/gdalrasterize.cpp
+++ b/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;
diff --git a/alg/gdalrasterpolygonenumerator.cpp b/alg/gdalrasterpolygonenumerator.cpp
index 4294d67c9cc8..a6fb42bca5a5 100644
--- a/alg/gdalrasterpolygonenumerator.cpp
+++ b/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>;
 
diff --git a/alg/gdalsievefilter.cpp b/alg/gdalsievefilter.cpp
index 6ac4368e8d03..540e20b22161 100644
--- a/alg/gdalsievefilter.cpp
+++ b/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      */
diff --git a/alg/gdalwarpkernel.cpp b/alg/gdalwarpkernel.cpp
index 503092e90a92..bc0dd28dbe77 100644
--- a/alg/gdalwarpkernel.cpp
+++ b/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]);
diff --git a/alg/polygonize.cpp b/alg/polygonize.cpp
index fc6f2e07d57b..8ba5c53344c7 100644
--- a/alg/polygonize.cpp
+++ b/alg/polygonize.cpp
@@ -406,14 +406,16 @@ static void AddEdges( GInt32 *panThisLineId, GInt32 *panLastLineId,
         if( nThisId != -1 )
         {
             if( papoPoly[nThisId] == nullptr )
-                papoPoly[nThisId] = new RPolygon( panPolyValue[nThisId] );
+                // FIXME loss of precision for [U]Int64
+                papoPoly[nThisId] = new RPolygon( static_cast<double>(panPolyValue[nThisId]) );
 
             papoPoly[nThisId]->AddSegment( iXReal, iY, iXReal+1, iY, 1 );
         }
         if( nPreviousId != -1 )
         {
             if( papoPoly[nPreviousId] == nullptr )
-                papoPoly[nPreviousId] = new RPolygon(panPolyValue[nPreviousId]);
+                // FIXME loss of precision for [U]Int64
+                papoPoly[nPreviousId] = new RPolygon( static_cast<double>(panPolyValue[nPreviousId]) );
 
             papoPoly[nPreviousId]->AddSegment( iXReal, iY, iXReal+1, iY, 0 );
         }
@@ -424,7 +426,8 @@ static void AddEdges( GInt32 *panThisLineId, GInt32 *panLastLineId,
         if( nThisId != -1 )
         {
             if( papoPoly[nThisId] == nullptr )
-                papoPoly[nThisId] = new RPolygon(panPolyValue[nThisId]);
+                // FIXME loss of precision for [U]Int64
+                papoPoly[nThisId] = new RPolygon( static_cast<double>(panPolyValue[nThisId]) );
 
             papoPoly[nThisId]->AddSegment( iXReal+1, iY, iXReal+1, iY+1, 1 );
         }
@@ -432,7 +435,8 @@ static void AddEdges( GInt32 *panThisLineId, GInt32 *panLastLineId,
         if( nRightId != -1 )
         {
             if( papoPoly[nRightId] == nullptr )
-                papoPoly[nRightId] = new RPolygon(panPolyValue[nRightId]);
+                // FIXME loss of precision for [U]Int64
+                papoPoly[nRightId] = new RPolygon( static_cast<double>(panPolyValue[nRightId]) );
 
             papoPoly[nRightId]->AddSegment( iXReal+1, iY, iXReal+1, iY+1, 0 );
         }
@@ -899,7 +903,7 @@ GBool GDALFloatEquals( float A, float B )
  * can be provided to determine which pixels are eligible for processing.
  *
  * Note that currently the source pixel band values are read into a
- * signed 32bit integer buffer (Int32), so floating point or complex
+ * signed 64bit integer buffer (Int64), so floating point or complex
  * bands will be implicitly truncated before processing. If you want to use a
  * version using 32bit float buffers, see GDALFPolygonize().
   *
@@ -954,14 +958,14 @@ GDALPolygonize( GDALRasterBandH hSrcBand,
                 void * pProgressArg )
 
 {
-    return GDALPolygonizeT<GInt32, IntEqualityTest>(hSrcBand,
+    return GDALPolygonizeT<std::int64_t, IntEqualityTest>(hSrcBand,
                                                     hMaskBand,
                                                     hOutLayer,
                                                     iPixValField,
                                                     papszOptions,
                                                     pfnProgress,
                                                     pProgressArg,
-                                                    GDT_Int32);
+                                                    GDT_Int64);
 }
 
 /************************************************************************/
diff --git a/apps/gdal_create.cpp b/apps/gdal_create.cpp
index 405785ace3d3..b55b78c3d49b 100644
--- a/apps/gdal_create.cpp
+++ b/apps/gdal_create.cpp
@@ -46,7 +46,7 @@ static void Usage()
             "       [-outsize xsize ysize]\n"
             "       [-bands count]\n"
             "       [-burn value]*\n"
-            "       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/\n"
+            "       [-ot {Byte/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/\n"
             "             CInt16/CInt32/CFloat32/CFloat64}] [-strict]\n"
             "       [-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]\n"
             "       [-mo \"META-TAG=VALUE\"]* [-q]\n"
diff --git a/apps/gdal_rasterize_bin.cpp b/apps/gdal_rasterize_bin.cpp
index 757849738a62..47351a52dbca 100644
--- a/apps/gdal_rasterize_bin.cpp
+++ b/apps/gdal_rasterize_bin.cpp
@@ -48,7 +48,7 @@ static void Usage(const char* pszErrorMsg = nullptr)
         "       [-dialect dialect] [-of format] [-a_srs srs_def] [-to \"NAME=VALUE\"]*\n"
         "       [-co \"NAME=VALUE\"]* [-a_nodata value] [-init value]*\n"
         "       [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]\n"
-        "       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/\n"
+        "       [-ot {Byte/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/\n"
         "             CInt16/CInt32/CFloat32/CFloat64}] [-optim {[AUTO]/VECTOR/RASTER}] [-q]\n"
         "       <src_datasource> <dst_filename>\n" );
 
diff --git a/apps/gdal_translate_bin.cpp b/apps/gdal_translate_bin.cpp
index 7868109350db..59b52d1170a1 100644
--- a/apps/gdal_translate_bin.cpp
+++ b/apps/gdal_translate_bin.cpp
@@ -46,7 +46,7 @@ static void Usage(const char* pszErrorMsg, int bShort)
 
 {
     printf( "Usage: gdal_translate [--help-general] [--long-usage]\n"
-            "       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/\n"
+            "       [-ot {Byte/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/\n"
             "             CInt16/CInt32/CFloat32/CFloat64}] [-strict]\n"
             "       [-if format]* [-of format]\n"
             "       [-b band] [-mask band] [-expand {gray|rgb|rgba}]\n"
diff --git a/apps/gdal_translate_lib.cpp b/apps/gdal_translate_lib.cpp
index 0381d96c33ba..2f4e9986aeb4 100644
--- a/apps/gdal_translate_lib.cpp
+++ b/apps/gdal_translate_lib.cpp
@@ -1783,39 +1783,53 @@ GDALDatasetH GDALTranslate( const char *pszDest, GDALDatasetH hSrcDataset,
                                 !GDALDataTypeIsComplex(eBandType));
                 if( bSrcIsInteger && bDstIsInteger )
                 {
-                    GInt32 nDstMin = 0;
-                    GUInt32 nDstMax = 0;
+                    std::int64_t nDstMin = 0;
+                    std::uint64_t nDstMax = 0;
                     switch( eBandType )
                     {
                         case GDT_Byte:
-                            nDstMin = 0;
-                            nDstMax = 255;
+                            nDstMin = std::numeric_limits<std::uint8_t>::min();
+                            nDstMax = std::numeric_limits<std::uint8_t>::max();
                             break;
                         case GDT_UInt16:
-                            nDstMin = 0;
-                            nDstMax = 65535;
+                            nDstMin = std::numeric_limits<std::uint16_t>::min();
+                            nDstMax = std::numeric_limits<std::uint16_t>::max();
                             break;
                         case GDT_Int16:
-                            nDstMin = -32768;
-                            nDstMax = 32767;
+                            nDstMin = std::numeric_limits<std::int16_t>::min();
+                            nDstMax = std::numeric_limits<std::int16_t>::max();
                             break;
                         case GDT_UInt32:
-                            nDstMin = 0;
-                            nDstMax = 0xFFFFFFFFU;
+                            nDstMin = std::numeric_limits<std::uint32_t>::min();
+                            nDstMax = std::numeric_limits<std::uint32_t>::max();
                             break;
                         case GDT_Int32:
-                            nDstMin = 0x80000000;
-                            nDstMax = 0x7FFFFFFF;
+                            nDstMin = std::numeric_limits<std::int32_t>::min();
+                            nDstMax = std::numeric_limits<std::int32_t>::max();
+                            break;
+                        case GDT_UInt64:
+                            nDstMin = std::numeric_limits<std::uint64_t>::min();
+                            nDstMax = std::numeric_limits<std::uint64_t>::max();
+                            break;
+                        case GDT_Int64:
+                            nDstMin = std::numeric_limits<std::int64_t>::min();
+                            nDstMax = std::numeric_limits<std::int64_t>::max();
                             break;
                         default:
                             CPLAssert(false);
                             break;
                     }
 
-                    GInt32 nMin = atoi(pszMin);
-                    GUInt32 nMax = static_cast<GUInt32>(strtoul(pszMax, nullptr, 10));
-                    if( nMin < nDstMin || nMax > nDstMax )
-                        bFilterOutStatsMetadata = true;
+                    try
+                    {
+                        const auto nMin = std::stoll(pszMin);
+                        const auto nMax = std::stoull(pszMax);
+                        if( nMin < nDstMin || nMax > nDstMax )
+                            bFilterOutStatsMetadata = true;
+                    }
+                    catch( const std::exception& )
+                    {
+                    }
                 }
                 // Float64 is large enough to hold all integer <= 32 bit or float32 values
                 // there might be other OK cases, but ere on safe side for now
diff --git a/autotest/cpp/test_cpl.cpp b/autotest/cpp/test_cpl.cpp
index c01324ab30fb..0b2c5c1b3eda 100644
--- a/autotest/cpp/test_cpl.cpp
+++ b/autotest/cpp/test_cpl.cpp
@@ -2714,12 +2714,12 @@ namespace tut
         }
         {
             CPLJSonStreamingWriter x(nullptr, nullptr);
-            x.Add(static_cast<GIntBig>(-10000) * 1000000);
+            x.Add(static_cast<std::int64_t>(-10000) * 1000000);
             ensure_equals( x.GetString(), std::string("-10000000000") );
         }
         {
             CPLJSonStreamingWriter x(nullptr, nullptr);
-            x.Add(static_cast<GUInt64>(10000) * 1000000);
+            x.Add(static_cast<std::uint64_t>(10000) * 1000000);
             ensure_equals( x.GetString(), std::string("10000000000") );
         }
         {
diff --git a/autotest/cpp/test_gdal.cpp b/autotest/cpp/test_gdal.cpp
index 7e2dc4e985ae..d1bfb015f22e 100644
--- a/autotest/cpp/test_gdal.cpp
+++ b/autotest/cpp/test_gdal.cpp
@@ -160,8 +160,13 @@ namespace tut
         }
 
         ENSURE_EQUALS(GDALDataTypeUnion(GDT_Int16, GDT_UInt16), GDT_Int32);
-        ENSURE_EQUALS(GDALDataTypeUnion(GDT_Int16, GDT_UInt32), GDT_Float64);
-        ENSURE_EQUALS(GDALDataTypeUnion(GDT_UInt32, GDT_Int16), GDT_Float64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_Int16, GDT_UInt32), GDT_Int64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_UInt32, GDT_Int16), GDT_Int64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_Int64, GDT_UInt64), GDT_Float64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_Int64, GDT_Float32), GDT_Float64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_Int64, GDT_Float64), GDT_Float64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_UInt64, GDT_Float32), GDT_Float64);
+        ENSURE_EQUALS(GDALDataTypeUnion(GDT_UInt64, GDT_Float64), GDT_Float64);
         ENSURE_EQUALS(GDALDataTypeUnion(GDT_UInt32, GDT_CInt16), GDT_CFloat64);
         ENSURE_EQUALS(GDALDataTypeUnion(GDT_Float32, GDT_CInt32), GDT_CFloat64);
         ENSURE_EQUALS(GDALDataTypeUnion(GDT_CInt16, GDT_UInt32), GDT_CFloat64);
@@ -204,6 +209,9 @@ namespace tut
         ENSURE_EQUALS(GDALFindDataType(64, false /* signed */, true /* floating */, false /* complex */), GDT_Float64);
         ENSURE_EQUALS(GDALFindDataType(64, false /* signed */, true /* floating */, true /* complex */), GDT_CFloat64);
 
+        ENSURE_EQUALS(GDALFindDataType(64, false /* signed */, false /* floating */, false /* complex */), GDT_UInt64);
+        ENSURE_EQUALS(GDALFindDataType(64, true /* signed */, false /* floating */, false /* complex */), GDT_Int64);
+
         ENSURE_EQUALS(GDALDataTypeUnionWithValue(GDT_Byte, -32768, 0), GDT_Int16);
         ENSURE_EQUALS(GDALDataTypeUnionWithValue(GDT_Byte, -32769, 0), GDT_Int32);
         ENSURE_EQUALS(GDALDataTypeUnionWithValue(GDT_Float32, -99999, 0), GDT_Float32);
@@ -243,6 +251,13 @@ namespace tut
         ensure( GDALAdjustValueToDataType(GDT_Int32,-10000000.0,&bClamped,&bRounded) == -10000000.0 && !bClamped && !bRounded);
         ensure( GDALAdjustValueToDataType(GDT_Int32,10000000.0,&bClamped,&bRounded) == 10000000.0 && !bClamped && !bRounded);
 
+        ensure( GDALAdjustValueToDataType(GDT_UInt64,10000000000.0,&bClamped,&bRounded) == 10000000000.0 && !bClamped && !bRounded);
+        ensure( GDALAdjustValueToDataType(GDT_UInt64,10000000000.4,&bClamped,&bRounded) == 10000000000.0 && !bClamped && bRounded);
+        ensure( GDALAdjustValueToDataType(GDT_UInt64,-1,&bClamped,&bRounded) == 0.0 && bClamped && !bRounded);
+
+        ensure( GDALAdjustValueToDataType(GDT_Int64,-10000000000.0,&bClamped,&bRounded) == -10000000000.0 && !bClamped && !bRounded);
+        ensure( GDALAdjustValueToDataType(GDT_Int64,10000000000.0,&bClamped,&bRounded) == 10000000000.0 && !bClamped && !bRounded);
+
         ensure( GDALAdjustValueToDataType(GDT_Float32, 0.0,&bClamped,&bRounded) == 0.0 && !bClamped && !bRounded);
         ensure( GDALAdjustValueToDataType(GDT_Float32, 1e-50,&bClamped,&bRounded) == 0.0 && !bClamped && !bRounded);
         ensure( GDALAdjustValueToDataType(GDT_Float32, 1.23,&bClamped,&bRounded) == static_cast<double>(1.23f) && !bClamped && !bRounded);
@@ -458,6 +473,8 @@ namespace tut
         ensure_equals( GDALDataTypeIsInteger(GDT_Int16), TRUE );
         ensure_equals( GDALDataTypeIsInteger(GDT_UInt32), TRUE );
         ensure_equals( GDALDataTypeIsInteger(GDT_Int32), TRUE );
+        ensure_equals( GDALDataTypeIsInteger(GDT_UInt64), TRUE );
+        ensure_equals( GDALDataTypeIsInteger(GDT_Int64), TRUE );
         ensure( !GDALDataTypeIsInteger(GDT_Float32) );
         ensure( !GDALDataTypeIsInteger(GDT_Float64) );
         ensure_equals( GDALDataTypeIsInteger(GDT_CInt16), TRUE );
@@ -475,6 +492,8 @@ namespace tut
         ensure( !GDALDataTypeIsFloating(GDT_Int16) );
         ensure( !GDALDataTypeIsFloating(GDT_UInt32) );
         ensure( !GDALDataTypeIsFloating(GDT_Int32) );
+        ensure( !GDALDataTypeIsFloating(GDT_UInt64) );
+        ensure( !GDALDataTypeIsFloating(GDT_Int64) );
         ensure_equals( GDALDataTypeIsFloating(GDT_Float32), TRUE );
         ensure_equals( GDALDataTypeIsFloating(GDT_Float64), TRUE );
         ensure( !GDALDataTypeIsFloating(GDT_CInt16) );
@@ -492,6 +511,8 @@ namespace tut
         ensure( !GDALDataTypeIsComplex(GDT_Int16) );
         ensure( !GDALDataTypeIsComplex(GDT_UInt32) );
         ensure( !GDALDataTypeIsComplex(GDT_Int32) );
+        ensure( !GDALDataTypeIsComplex(GDT_UInt64) );
+        ensure( !GDALDataTypeIsComplex(GDT_Int64) );
         ensure( !GDALDataTypeIsComplex(GDT_Float32) );
         ensure( !GDALDataTypeIsComplex(GDT_Float64) );
         ensure_equals( GDALDataTypeIsComplex(GDT_CInt16), TRUE );
@@ -508,6 +529,8 @@ namespace tut
         ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_Int16) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_UInt32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_Int32) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_UInt64) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_Int64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_Float64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Byte, GDT_CInt16) );
@@ -520,6 +543,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_Int16) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_UInt32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_Int32) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_UInt64) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_Int64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_UInt16, GDT_CInt16) );
@@ -532,6 +557,8 @@ namespace tut
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int16, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Int16, GDT_UInt32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int16, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int16, GDT_UInt64) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_Int16, GDT_Int64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int16, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int16, GDT_Float64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int16, GDT_CInt16) );
@@ -544,6 +571,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_Int16) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_Int32) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_UInt64) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_UInt32, GDT_CInt16) );
@@ -556,6 +585,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_Int32, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Int32, GDT_UInt32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int32, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int32, GDT_UInt64) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_Int32, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Int32, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int32, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Int32, GDT_CInt16) );
@@ -563,11 +594,41 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_Int32, GDT_CFloat32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Int32, GDT_CFloat64) );
 
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_Byte) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_UInt16) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_Int16) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_UInt32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_Int32) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_Int64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_Float32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_Float64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_CInt16) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_CInt32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_CFloat32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_UInt64, GDT_CFloat64) );
+
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_Byte) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_UInt16) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_Int16) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_UInt32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_UInt64) );
+        ensure( !GDALDataTypeIsConversionLossy(GDT_Int64, GDT_Int64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_Float32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_Float64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_CInt16) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_CInt32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_CFloat32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Int64, GDT_CFloat64) );
+
         ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_Byte) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_UInt16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_Int64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Float32, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Float32, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float32, GDT_CInt16) );
@@ -580,6 +641,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_Float32) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_Float64, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_Float64, GDT_CInt16) );
@@ -592,6 +655,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_Float32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_Float64) );
         ensure( !GDALDataTypeIsConversionLossy(GDT_CInt16, GDT_CInt16) );
@@ -604,6 +669,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_Float32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CInt32, GDT_CInt16) );
@@ -616,6 +683,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_Float32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat32, GDT_CInt16) );
@@ -628,6 +697,8 @@ namespace tut
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_Int16) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_UInt32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_Int32) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_UInt64) );
+        ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_Int64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_Float32) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_Float64) );
         ensure( GDALDataTypeIsConversionLossy(GDT_CFloat64, GDT_CInt16) );
diff --git a/autotest/cpp/testcopywords.cpp b/autotest/cpp/testcopywords.cpp
index 497f0aae5fd2..6ba6cd5133f1 100644
--- a/autotest/cpp/testcopywords.cpp
+++ b/autotest/cpp/testcopywords.cpp
@@ -125,6 +125,10 @@ template <class InType, class ConstantType> void FromR_2(GDALDataType intype, Co
         Test<InType,GInt32,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
     else if (outtype == GDT_UInt32)
         Test<InType,GUInt32,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
+    else if (outtype == GDT_Int64)
+        Test<InType,std::int64_t,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
+    else if (outtype == GDT_UInt64)
+        Test<InType,std::uint64_t,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
     else if (outtype == GDT_Float32)
         Test<InType,float,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
     else if (outtype == GDT_Float64)
@@ -152,6 +156,10 @@ void FromR(GDALDataType intype, ConstantType inval, ConstantType invali, GDALDat
         FromR_2<GInt32,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
     else if (intype == GDT_UInt32)
         FromR_2<GUInt32,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
+    else if (intype == GDT_Int64)
+        FromR_2<std::int64_t,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
+    else if (intype == GDT_UInt64)
+        FromR_2<std::uint64_t,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
     else if (intype == GDT_Float32)
         FromR_2<float,ConstantType>(intype, inval, invali, outtype, outval, outvali, numLine);
     else if (intype == GDT_Float64)
@@ -173,7 +181,7 @@ void FromR(GDALDataType intype, ConstantType inval, ConstantType invali, GDALDat
 #define FROM_C(intype, inval, invali, outtype, outval, outvali) FromR<GIntBig>(intype, inval, invali, outtype, outval, outvali, __LINE__)
 #define FROM_C_F(intype, inval, invali, outtype, outval, outvali) FromR<double>(intype, inval, invali, outtype, outval, outvali, __LINE__)
 
-#define IS_UNSIGNED(x) (x == GDT_Byte || x == GDT_UInt16 || x == GDT_UInt32)
+#define IS_UNSIGNED(x) (x == GDT_Byte || x == GDT_UInt16 || x == GDT_UInt32 || x == GDT_UInt64)
 #define IS_FLOAT(x) (x == GDT_Float32 || x == GDT_Float64 || x == GDT_CFloat32 || x == GDT_CFloat64)
 
 #define CST_3000000000 (((GIntBig)3000) * 1000 * 1000)
@@ -182,7 +190,7 @@ void FromR(GDALDataType intype, ConstantType inval, ConstantType invali, GDALDat
 static void check_GDT_Byte()
 {
     /* GDT_Byte */
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype<GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_R(GDT_Byte, 0, outtype, 0);
         FROM_R(GDT_Byte, 127, outtype, 127);
@@ -222,13 +230,15 @@ static void check_GDT_Int16()
     FROM_R(GDT_Int16, -32000, GDT_UInt16, 0); /* clamp */
     FROM_R(GDT_Int16, -32000, GDT_Int32, -32000);
     FROM_R(GDT_Int16, -32000, GDT_UInt32, 0); /* clamp */
+    FROM_R(GDT_Int16, -32000, GDT_Int64, -32000);
+    FROM_R(GDT_Int16, -32000, GDT_UInt64, 0); /* clamp */
     FROM_R(GDT_Int16, -32000, GDT_Float32, -32000);
     FROM_R(GDT_Int16, -32000, GDT_Float64, -32000);
     FROM_R(GDT_Int16, -32000, GDT_CInt16, -32000);
     FROM_R(GDT_Int16, -32000, GDT_CInt32, -32000);
     FROM_R(GDT_Int16, -32000, GDT_CFloat32, -32000);
     FROM_R(GDT_Int16, -32000, GDT_CFloat64, -32000);
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_R(GDT_Int16, 127, outtype, 127);
     }
@@ -238,6 +248,8 @@ static void check_GDT_Int16()
     FROM_R(GDT_Int16, 32000, GDT_UInt16, 32000);
     FROM_R(GDT_Int16, 32000, GDT_Int32, 32000);
     FROM_R(GDT_Int16, 32000, GDT_UInt32, 32000);
+    FROM_R(GDT_Int16, 32000, GDT_Int64, 32000);
+    FROM_R(GDT_Int16, 32000, GDT_UInt64, 32000);
     FROM_R(GDT_Int16, 32000, GDT_Float32, 32000);
     FROM_R(GDT_Int16, 32000, GDT_Float64, 32000);
     FROM_R(GDT_Int16, 32000, GDT_CInt16, 32000);
@@ -249,7 +261,7 @@ static void check_GDT_Int16()
 static void check_GDT_UInt16()
 {
     /* GDT_UInt16 */
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_R(GDT_UInt16, 0, outtype, 0);
         FROM_R(GDT_UInt16, 127, outtype, 127);
@@ -260,6 +272,8 @@ static void check_GDT_UInt16()
     FROM_R(GDT_UInt16, 65000, GDT_UInt16, 65000);
     FROM_R(GDT_UInt16, 65000, GDT_Int32, 65000);
     FROM_R(GDT_UInt16, 65000, GDT_UInt32, 65000);
+    FROM_R(GDT_UInt16, 65000, GDT_Int64, 65000);
+    FROM_R(GDT_UInt16, 65000, GDT_UInt64, 65000);
     FROM_R(GDT_UInt16, 65000, GDT_Float32, 65000);
     FROM_R(GDT_UInt16, 65000, GDT_Float64, 65000);
     FROM_R(GDT_UInt16, 65000, GDT_CInt16, 32767); /* clamp */
@@ -274,15 +288,17 @@ static void check_GDT_Int32()
     FROM_R(GDT_Int32, -33000, GDT_Byte, 0); /* clamp */
     FROM_R(GDT_Int32, -33000, GDT_Int16, -32768); /* clamp */
     FROM_R(GDT_Int32, -33000, GDT_UInt16, 0); /* clamp */
-    FROM_R(GDT_Int32, -33000, GDT_Int32, -33000); /* clamp */
+    FROM_R(GDT_Int32, -33000, GDT_Int32, -33000);
     FROM_R(GDT_Int32, -33000, GDT_UInt32, 0); /* clamp */
+    FROM_R(GDT_Int32, -33000, GDT_Int64, -33000);
+    FROM_R(GDT_Int32, -33000, GDT_UInt64, 0); /* clamp */
     FROM_R(GDT_Int32, -33000, GDT_Float32, -33000);
     FROM_R(GDT_Int32, -33000, GDT_Float64, -33000);
     FROM_R(GDT_Int32, -33000, GDT_CInt16, -32768); /* clamp */
     FROM_R(GDT_Int32, -33000, GDT_CInt32, -33000);
     FROM_R(GDT_Int32, -33000, GDT_CFloat32, -33000);
     FROM_R(GDT_Int32, -33000, GDT_CFloat64, -33000);
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_R(GDT_Int32, 127, outtype, 127);
     }
@@ -292,6 +308,8 @@ static void check_GDT_Int32()
     FROM_R(GDT_Int32, 67000, GDT_UInt16, 65535);  /* clamp */
     FROM_R(GDT_Int32, 67000, GDT_Int32, 67000);
     FROM_R(GDT_Int32, 67000, GDT_UInt32, 67000);
+    FROM_R(GDT_Int32, 67000, GDT_Int64, 67000);
+    FROM_R(GDT_Int32, 67000, GDT_UInt64, 67000);
     FROM_R(GDT_Int32, 67000, GDT_Float32, 67000);
     FROM_R(GDT_Int32, 67000, GDT_Float64, 67000);
     FROM_R(GDT_Int32, 67000, GDT_CInt16, 32767);  /* clamp */
@@ -303,7 +321,7 @@ static void check_GDT_Int32()
 static void check_GDT_UInt32()
 {
     /* GDT_UInt32 */
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_R(GDT_UInt32, 0, outtype, 0);
         FROM_R(GDT_UInt32, 127, outtype, 127);
@@ -314,6 +332,8 @@ static void check_GDT_UInt32()
     FROM_R(GDT_UInt32, 3000000000U, GDT_UInt16, 65535);  /* clamp */
     FROM_R(GDT_UInt32, 3000000000U, GDT_Int32, 2147483647);  /* clamp */
     FROM_R(GDT_UInt32, 3000000000U, GDT_UInt32, 3000000000U);
+    FROM_R(GDT_UInt32, 3000000000U, GDT_Int64, 3000000000U);
+    FROM_R(GDT_UInt32, 3000000000U, GDT_UInt64, 3000000000U);
     FROM_R(GDT_UInt32, 3000000000U, GDT_Float32, 3000000000U);
     FROM_R(GDT_UInt32, 3000000000U, GDT_Float64, 3000000000U);
     FROM_R(GDT_UInt32, 3000000000U, GDT_CInt16, 32767);  /* clamp */
@@ -322,13 +342,74 @@ static void check_GDT_UInt32()
     FROM_R(GDT_UInt32, 3000000000U, GDT_CFloat64, 3000000000U);
 }
 
+static void check_GDT_Int64()
+{
+    /* GDT_Int64 */
+    FROM_R(GDT_Int64, -33000, GDT_Byte, 0); /* clamp */
+    FROM_R(GDT_Int64, -33000, GDT_Int16, -32768); /* clamp */
+    FROM_R(GDT_Int64, -33000, GDT_UInt16, 0); /* clamp */
+    FROM_R(GDT_Int64, -33000, GDT_Int32, -33000);
+    FROM_R(GDT_Int64, -33000, GDT_UInt32, 0); /* clamp */
+    FROM_R(GDT_Int64, -33000, GDT_Int64, -33000);
+    FROM_R(GDT_Int64, -33000, GDT_UInt64, 0); /* clamp */
+    FROM_R(GDT_Int64, -33000, GDT_Float32, -33000);
+    FROM_R(GDT_Int64, -33000, GDT_Float64, -33000);
+    FROM_R(GDT_Int64, -33000, GDT_CInt16, -32768); /* clamp */
+    FROM_R(GDT_Int64, -33000, GDT_CInt32, -33000);
+    FROM_R(GDT_Int64, -33000, GDT_CFloat32, -33000);
+    FROM_R(GDT_Int64, -33000, GDT_CFloat64, -33000);
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
+    {
+        FROM_R(GDT_Int64, 127, outtype, 127);
+    }
+
+    FROM_R(GDT_Int64, 67000, GDT_Byte, 255); /* clamp */
+    FROM_R(GDT_Int64, 67000, GDT_Int16, 32767);  /* clamp */
+    FROM_R(GDT_Int32, 67000, GDT_UInt16, 65535);  /* clamp */
+    FROM_R(GDT_Int32, 67000, GDT_Int32, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_UInt32, 67000);
+    FROM_R(GDT_Int32, 67000, GDT_Int64, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_UInt64, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_Float32, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_Float64, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_CInt16, 32767);  /* clamp */
+    FROM_R(GDT_Int64, 67000, GDT_CInt32, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_CFloat32, 67000);
+    FROM_R(GDT_Int64, 67000, GDT_CFloat64, 67000);
+}
+
+static void check_GDT_UInt64()
+{
+    /* GDT_UInt64 */
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
+    {
+        FROM_R(GDT_UInt64, 0, outtype, 0);
+        FROM_R(GDT_UInt64, 127, outtype, 127);
+    }
+
+    std::uint64_t nVal = static_cast<std::uint64_t>(3000000000) * 1000;
+    FROM_R(GDT_UInt64, nVal, GDT_Byte, 255); /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_Int16, 32767);  /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_UInt16, 65535);  /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_Int32, 2147483647);  /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_UInt32, 4294967295U); /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_Int64, nVal);
+    FROM_R(GDT_UInt64, nVal, GDT_UInt64, nVal);
+    FROM_R(GDT_UInt64, nVal, GDT_Float32, static_cast<uint64_t>(static_cast<float>(nVal)));
+    FROM_R(GDT_UInt64, nVal, GDT_Float64, nVal);
+    FROM_R(GDT_UInt64, nVal, GDT_CInt16, 32767);  /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_CInt32, 2147483647);  /* clamp */
+    FROM_R(GDT_UInt64, nVal, GDT_CFloat32, static_cast<uint64_t>(static_cast<float>(nVal)));
+    FROM_R(GDT_UInt64, nVal, GDT_CFloat64, nVal);
+}
+
 static void check_GDT_Float32and64()
 {
     /* GDT_Float32 and GDT_Float64 */
     for(int i=0;i<2;i++)
     {
         GDALDataType intype = (i == 0) ? GDT_Float32 : GDT_Float64;
-        for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+        for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
         {
             if (IS_FLOAT(outtype))
             {
@@ -404,7 +485,7 @@ static void check_GDT_CInt16()
     FROM_C(GDT_CInt16, -32000, -32500, GDT_CInt32, -32000, -32500);
     FROM_C(GDT_CInt16, -32000, -32500, GDT_CFloat32, -32000, -32500);
     FROM_C(GDT_CInt16, -32000, -32500, GDT_CFloat64, -32000, -32500);
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_C(GDT_CInt16, 127, 128, outtype, 127, 128);
     }
@@ -436,7 +517,7 @@ static void check_GDT_CInt32()
     FROM_C(GDT_CInt32, -33000, -33500, GDT_CInt32, -33000, -33500);
     FROM_C(GDT_CInt32, -33000, -33500, GDT_CFloat32, -33000, -33500);
     FROM_C(GDT_CInt32, -33000, -33500, GDT_CFloat64, -33000, -33500);
-    for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+    for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
     {
         FROM_C(GDT_CInt32, 127, 128, outtype, 127, 128);
     }
@@ -460,7 +541,7 @@ static void check_GDT_CFloat32and64()
     for(int i=0;i<2;i++)
     {
         GDALDataType intype = (i == 0) ? GDT_CFloat32 : GDT_CFloat64;
-        for(GDALDataType outtype=GDT_Byte; outtype<=GDT_CFloat64;outtype = (GDALDataType)(outtype + 1))
+        for(GDALDataType outtype=GDT_Byte; outtype< GDT_TypeCount;outtype = (GDALDataType)(outtype + 1))
         {
             if (IS_FLOAT(outtype))
             {
@@ -579,6 +660,8 @@ void CheckPacked(GDALDataType eIn, GDALDataType eOut)
         case GDT_Int16: CheckPacked<Tin, GInt16>(eIn, eOut); break;
         case GDT_UInt32: CheckPacked<Tin, GUInt32>(eIn, eOut); break;
         case GDT_Int32: CheckPacked<Tin, GInt32>(eIn, eOut); break;
+        case GDT_UInt64: CheckPacked<Tin, std::uint64_t>(eIn, eOut); break;
+        case GDT_Int64: CheckPacked<Tin, std::int64_t>(eIn, eOut); break;
         case GDT_Float32: CheckPacked<Tin, float>(eIn, eOut); break;
         case GDT_Float64: CheckPacked<Tin, double>(eIn, eOut); break;
         default:
@@ -597,6 +680,8 @@ static void CheckPacked(GDALDataType eIn, GDALDataType eOut)
         case GDT_Int16: CheckPacked<GInt16>(eIn, eOut); break;
         case GDT_UInt32: CheckPacked<GUInt32>(eIn, eOut); break;
         case GDT_Int32: CheckPacked<GInt32>(eIn, eOut); break;
+        case GDT_UInt64: CheckPacked<std::uint64_t>(eIn, eOut); break;
+        case GDT_Int64: CheckPacked<std::int64_t>(eIn, eOut); break;
         case GDT_Float32: CheckPacked<float>(eIn, eOut); break;
         case GDT_Float64: CheckPacked<double>(eIn, eOut); break;
         default:
@@ -615,6 +700,8 @@ int main(int /* argc */, char* /* argv */ [])
     check_GDT_UInt16();
     check_GDT_Int32();
     check_GDT_UInt32();
+    check_GDT_Int64();
+    check_GDT_UInt64();
     check_GDT_Float32and64();
     check_GDT_CInt16();
     check_GDT_CInt32();
@@ -704,10 +791,12 @@ int main(int /* argc */, char* /* argv */ [])
     free(pIn);
     free(pOut);
 
-    for( GDALDataType eIn = GDT_Byte; eIn <= GDT_Float64; eIn = static_cast<GDALDataType>(eIn + 1) )
+    for( GDALDataType eIn = GDT_Byte; eIn < GDT_TypeCount; eIn = static_cast<GDALDataType>(eIn + 1) )
     {
-        for( GDALDataType eOut = GDT_Byte; eOut <= GDT_Float64; eOut = static_cast<GDALDataType>(eOut + 1) )
+        if( GDALDataTypeIsComplex(eIn) ) continue;
+        for( GDALDataType eOut = GDT_Byte; eOut < GDT_TypeCount; eOut = static_cast<GDALDataType>(eOut + 1) )
         {
+            if( GDALDataTypeIsComplex(eOut) ) continue;
             CheckPacked(eIn, eOut);
         }
     }
diff --git a/autotest/cpp/testperfcopywords.cpp b/autotest/cpp/testperfcopywords.cpp
index b564be4a6a0d..431bb2109156 100644
--- a/autotest/cpp/testperfcopywords.cpp
+++ b/autotest/cpp/testperfcopywords.cpp
@@ -44,9 +44,9 @@ int main(int /* argc */, char* /* argv */ [])
 
     clock_t start, end;
 
-    for(intype=GDT_Byte; intype<=GDT_CFloat64;intype++)
+    for(intype=GDT_Byte; intype<GDT_TypeCount;intype++)
     {
-        for(outtype=GDT_Byte;outtype<=GDT_CFloat64;outtype++)
+        for(outtype=GDT_Byte;outtype<GDT_TypeCount;outtype++)
         {
             start = clock();
 
diff --git a/autotest/gcore/data/gtiff/int64.tif b/autotest/gcore/data/gtiff/int64.tif
new file mode 100644
index 000000000000..bac5e5e75b6a
Binary files /dev/null and b/autotest/gcore/data/gtiff/int64.tif differ
diff --git a/autotest/gcore/data/gtiff/uint64.tif b/autotest/gcore/data/gtiff/uint64.tif
new file mode 100644
index 000000000000..6eb73d77f09c
Binary files /dev/null and b/autotest/gcore/data/gtiff/uint64.tif differ
diff --git a/autotest/gcore/numpy_rw.py b/autotest/gcore/numpy_rw.py
index 0d6dc1e5d7aa..c314377fdb12 100755
--- a/autotest/gcore/numpy_rw.py
+++ b/autotest/gcore/numpy_rw.py
@@ -248,6 +248,8 @@ def test_numpy_rw_11():
                    ('int16', gdal.GDT_Int16, numpy.int16, -32767),
                    ('uint32', gdal.GDT_UInt32, numpy.uint32, 4294967295),
                    ('int32', gdal.GDT_Int32, numpy.int32, -2147483648),
+                   ('uint64', gdal.GDT_UInt64, numpy.uint64, 4294967295 * 1000),
+                   ('int64', gdal.GDT_Int64, numpy.int64, -2147483648 * 1000),
                    ('float32', gdal.GDT_Float32, numpy.float32, 1.23),
                    ('float64', gdal.GDT_Float64, numpy.float64, 1.23456789),
                    ('cint16', gdal.GDT_CInt16, numpy.complex64, -32768 + 32767j),
@@ -277,6 +279,7 @@ def test_numpy_rw_11():
         ds = None
 
         ds = gdal.Open('/vsimem/' + type_tuple[0])
+        assert ds.GetRasterBand(1).DataType == type_tuple[1]
         ar2 = ds.ReadAsArray()
         ar3 = numpy.empty_like(ar2)
         ds.GetRasterBand(1).ReadAsArray(buf_obj=ar3)
@@ -331,7 +334,7 @@ def test_numpy_rw_13():
     ds.GetRasterBand(1).WriteArray(ar)
 
     # Try reading into unsupported array type
-    ar = numpy.empty([1, 2], dtype=numpy.int64)
+    ar = numpy.empty([1, 2], dtype=numpy.bool_)
     with pytest.raises(Exception, match='array does not have '
                              'corresponding GDAL data type'):
         ds.GetRasterBand(1).ReadAsArray(buf_obj=ar)
@@ -403,7 +406,7 @@ def test_numpy_rw_13():
     for i in range(3):
         ds.GetRasterBand(i + 1).WriteArray(ar[i])
 
-    ar = numpy.empty([3, 1, 2], dtype=numpy.int64)
+    ar = numpy.empty([3, 1, 2], dtype=numpy.bool_)
     with pytest.raises(Exception, match='array does not have '
                              'corresponding GDAL data type'):
         ds.ReadAsArray(buf_obj=ar)
diff --git a/autotest/gcore/tiff_ovr.py b/autotest/gcore/tiff_ovr.py
index ed24e2dd96df..c6dd800c4a19 100755
--- a/autotest/gcore/tiff_ovr.py
+++ b/autotest/gcore/tiff_ovr.py
@@ -2184,6 +2184,47 @@ def test_tiff_ovr_fallback_to_multiband_overview_generate():
 
     gdal.GetDriverByName('GTiff').Delete(filename)
 
+
+###############################################################################
+
+
+def test_tiff_ovr_int64():
+
+    temp_path = '/vsimem/test.tif'
+    ds = gdal.GetDriverByName('GTiff').Create(temp_path, 2, 1, 1, gdal.GDT_Int64)
+    ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1,
+                                    struct.pack('q' * 2, -10000000000, -10000000000))
+    del ds
+    ds = gdal.OpenEx(temp_path, gdal.GA_ReadOnly)
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_Int64
+    assert ds.BuildOverviews('nearest', overviewlist=[2]) == 0
+    del ds
+    ds = gdal.OpenEx(temp_path, gdal.GA_ReadOnly)
+    assert struct.unpack('q', ds.GetRasterBand(1).GetOverview(0).ReadRaster()) == (-10000000000,)
+    del ds
+    gdal.GetDriverByName('GTiff').Delete(temp_path)
+
+
+###############################################################################
+
+
+def test_tiff_ovr_uint64():
+
+    temp_path = '/vsimem/test.tif'
+    ds = gdal.GetDriverByName('GTiff').Create(temp_path, 2, 1, 1, gdal.GDT_UInt64)
+    ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1,
+                                    struct.pack('Q' * 2, 10000000000, 10000000000))
+    del ds
+    ds = gdal.OpenEx(temp_path, gdal.GA_ReadOnly)
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_UInt64
+    assert ds.BuildOverviews('nearest', overviewlist=[2]) == 0
+    del ds
+    ds = gdal.OpenEx(temp_path, gdal.GA_ReadOnly)
+    assert struct.unpack('Q', ds.GetRasterBand(1).GetOverview(0).ReadRaster()) == (10000000000,)
+    del ds
+    gdal.GetDriverByName('GTiff').Delete(temp_path)
+
+
 ###############################################################################
 # Cleanup
 
diff --git a/autotest/gcore/tiff_read.py b/autotest/gcore/tiff_read.py
index bc58836da2c3..388efd67ee6c 100755
--- a/autotest/gcore/tiff_read.py
+++ b/autotest/gcore/tiff_read.py
@@ -67,7 +67,9 @@
     ('seperate_strip.tif', 2, 15234),  # TODO: Spelling.
     ('contig_tiled.tif', 2, 15234),
     ('contig_strip.tif', 2, 15234),
-    ('empty1bit.tif', 1, 0)
+    ('empty1bit.tif', 1, 0),
+    ('gtiff/int64.tif', 1, 65535),
+    ('gtiff/uint64.tif', 1, 1),
 ]
 
 
diff --git a/autotest/gcore/tiff_write.py b/autotest/gcore/tiff_write.py
index b24b3fc013ba..8c85eb384d25 100755
--- a/autotest/gcore/tiff_write.py
+++ b/autotest/gcore/tiff_write.py
@@ -8312,6 +8312,22 @@ def test_tiff_write_predictor_2_float64():
     ds = None
     gdal.Unlink(filename)
 
+###############################################################################
+
+
+def test_tiff_write_uint64():
+
+    ut = gdaltest.GDALTest('GTiff', 'gtiff/uint64.tif', 1, 1)
+    return ut.testCreateCopy()
+
+###############################################################################
+
+
+def test_tiff_write_int64():
+
+    ut = gdaltest.GDALTest('GTiff', 'gtiff/int64.tif', 1, 65535)
+    return ut.testCreateCopy()
+
 
 def test_tiff_write_cleanup():
     gdaltest.tiff_drv = None
diff --git a/autotest/gdrivers/data/netcdf/int64.nc b/autotest/gdrivers/data/netcdf/int64.nc
new file mode 100644
index 000000000000..d97e37df15e2
Binary files /dev/null and b/autotest/gdrivers/data/netcdf/int64.nc differ
diff --git a/autotest/gdrivers/data/netcdf/uint64.nc b/autotest/gdrivers/data/netcdf/uint64.nc
new file mode 100644
index 000000000000..50412d873593
Binary files /dev/null and b/autotest/gdrivers/data/netcdf/uint64.nc differ
diff --git a/autotest/gdrivers/hdf5multidim.py b/autotest/gdrivers/hdf5multidim.py
index be8b4b9f2d84..51fff93e06c8 100755
--- a/autotest/gdrivers/hdf5multidim.py
+++ b/autotest/gdrivers/hdf5multidim.py
@@ -110,8 +110,8 @@ def test_hdf5_multidim_var_alldatatypes():
                       ('short_var', gdal.GDT_Int16, (-32768, -32767)),
                       ('uint_var', gdal.GDT_UInt32, (4294967294, 4294967293)),
                       ('int_var', gdal.GDT_Int32, (-2147483648, -2147483647)),
-                      ('uint64_var', gdal.GDT_Float64, (1.8446744073709552e+19, 1.8446744073709552e+19)),
-                      ('int64_var', gdal.GDT_Float64, (-9.223372036854776e+18, -9.223372036854776e+18)),
+                      ('uint64_var', gdal.GDT_UInt64, (18446744073709551613, 18446744073709551612)),
+                      ('int64_var', gdal.GDT_Int64, (-9223372036854775808, -9223372036854775807)),
                       ('float_var', gdal.GDT_Float32, (1.25, 2.25)),
                       ('double_var', gdal.GDT_Float64, (1.25125, 2.25125)),
                       ('complex_int16_var', gdal.GDT_CInt16, (-32768, -32767, -32766, -32765)),
@@ -134,6 +134,10 @@ def test_hdf5_multidim_var_alldatatypes():
             assert struct.unpack('I' * len(val), var.Read()) == val
         if dt == gdal.GDT_Int32:
             assert struct.unpack('i' * len(val), var.Read()) == val
+        if dt == gdal.GDT_UInt64:
+            assert struct.unpack('Q' * len(val), var.Read()) == val
+        if dt == gdal.GDT_Int64:
+            assert struct.unpack('q' * len(val), var.Read()) == val
         if dt == gdal.GDT_Float32:
             assert struct.unpack('f' * len(val), var.Read()) == val
         if dt == gdal.GDT_Float64:
diff --git a/autotest/gdrivers/netcdf.py b/autotest/gdrivers/netcdf.py
index 45eac6a197b6..c43a468aacfa 100755
--- a/autotest/gdrivers/netcdf.py
+++ b/autotest/gdrivers/netcdf.py
@@ -4782,9 +4782,6 @@ def test_netcdf_open_coords_no_georef_indexing_variables():
 
 def test_netcdf_metadata_sentinel5():
 
-    if not gdaltest.netcdf_drv_has_nc4:
-        pytest.skip()
-
     ds = gdal.Open('data/netcdf/fake_ISO_METADATA.nc')
     assert ds is not None
     assert "json:ISO_METADATA" in ds.GetMetadataDomainList()
@@ -5017,6 +5014,68 @@ def test_netcdf_read_gmt_file():
     assert gt == pytest.approx((-34.6671666666667, 0.001, 0.0, 35.58483333333329, 0.0, -0.001)), gt
 
 
+###############################################################################
+
+def test_netcdf_read_int64():
+
+    if not gdaltest.netcdf_drv_has_nc4:
+        pytest.skip()
+
+    ds = gdal.Open('data/netcdf/int64.nc')
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_Int64
+    assert struct.unpack('q' * 4, ds.ReadRaster()) == (10000000001, 1,
+                                                       -10000000000, 10000000000)
+
+
+###############################################################################
+
+
+def test_netcdf_write_int64():
+
+    if not gdaltest.netcdf_drv_has_nc4:
+        pytest.skip()
+
+    src_ds = gdal.Open('data/netcdf/int64.nc')
+    gdaltest.netcdf_drv.CreateCopy('tmp/int64.nc', src_ds)
+    ds = gdal.Open('tmp/int64.nc')
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_Int64
+    assert struct.unpack('q' * 4, ds.ReadRaster()) == (10000000001, 1,
+                                                       -10000000000, 10000000000)
+    ds = None
+    os.unlink('tmp/int64.nc')
+
+###############################################################################
+
+
+def test_netcdf_read_uint64():
+
+    if not gdaltest.netcdf_drv_has_nc4:
+        pytest.skip()
+
+    ds = gdal.Open('data/netcdf/uint64.nc')
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_UInt64
+    assert struct.unpack('Q' * 4, ds.ReadRaster()) == (10000000001, 1,
+                                                       0, 10000000000)
+
+
+###############################################################################
+
+
+def test_netcdf_write_uint64():
+
+    if not gdaltest.netcdf_drv_has_nc4:
+        pytest.skip()
+
+    src_ds = gdal.Open('data/netcdf/uint64.nc')
+    gdaltest.netcdf_drv.CreateCopy('tmp/uint64.nc', src_ds)
+    ds = gdal.Open('tmp/uint64.nc')
+    assert ds.GetRasterBand(1).DataType == gdal.GDT_UInt64
+    assert struct.unpack('Q' * 4, ds.ReadRaster()) == (10000000001, 1,
+                                                       0, 10000000000)
+    ds = None
+    os.unlink('tmp/uint64.nc')
+
+
 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/autotest/gdrivers/netcdf_multidim.py b/autotest/gdrivers/netcdf_multidim.py
index 878b9e6298c6..04b67b15c760 100755
--- a/autotest/gdrivers/netcdf_multidim.py
+++ b/autotest/gdrivers/netcdf_multidim.py
@@ -199,8 +199,8 @@ def test_netcdf_multidim_var_alldatatypes():
                       ('short_var', gdal.GDT_Int16, (-32768, -32767)),
                       ('uint_var', gdal.GDT_UInt32, (4294967294, 4294967293)),
                       ('int_var', gdal.GDT_Int32, (-2147483648, -2147483647)),
-                      ('uint64_var', gdal.GDT_Float64, (1.8446744073709552e+19, 1.8446744073709552e+19)),
-                      ('int64_var', gdal.GDT_Float64, (-9.223372036854776e+18, -9.223372036854776e+18)),
+                      ('uint64_var', gdal.GDT_UInt64, (18446744073709551613, 18446744073709551612)),
+                      ('int64_var', gdal.GDT_Int64, (-9223372036854775808, -9223372036854775807)),
                       ('float_var', gdal.GDT_Float32, (1.25, 2.25)),
                       ('double_var', gdal.GDT_Float64, (1.25125, 2.25125)),
                       ('complex_int16_var', gdal.GDT_CInt16, (-32768, -32767, -32766, -32765)),
@@ -223,6 +223,10 @@ def test_netcdf_multidim_var_alldatatypes():
             assert struct.unpack('I' * len(val), var.Read()) == val
         if dt == gdal.GDT_Int32:
             assert struct.unpack('i' * len(val), var.Read()) == val
+        if dt == gdal.GDT_UInt64:
+            assert struct.unpack('Q' * len(val), var.Read()) == val
+        if dt == gdal.GDT_Int64:
+            assert struct.unpack('q' * len(val), var.Read()) == val
         if dt == gdal.GDT_Float32:
             assert struct.unpack('f' * len(val), var.Read()) == val
         if dt == gdal.GDT_Float64:
@@ -799,8 +803,11 @@ def dims_from_non_netcdf(rg):
 
         dims_from_non_netcdf(rg)
 
-        for dt in (gdal.GDT_Byte, gdal.GDT_Int16, gdal.GDT_UInt16,
-                   gdal.GDT_Int32, gdal.GDT_UInt32, gdal.GDT_Float32,
+        for dt in (gdal.GDT_Byte,
+                   gdal.GDT_Int16, gdal.GDT_UInt16,
+                   gdal.GDT_Int32, gdal.GDT_UInt32,
+                   gdal.GDT_Int64, gdal.GDT_UInt64,
+                   gdal.GDT_Float32,
                    gdal.GDT_Float64, gdal.GDT_CInt16, gdal.GDT_CInt32,
                    gdal.GDT_CFloat32, gdal.GDT_CFloat64):
 
@@ -839,15 +846,15 @@ def dims_from_non_netcdf(rg):
 
         var = rg.CreateMDArray('var_as_nc_int64', [],
                 gdal.ExtendedDataType.Create(gdal.GDT_Float64), ['NC_TYPE=NC_INT64'])
-        assert var.GetDataType().GetNumericDataType() == gdal.GDT_Float64
-        assert var.Write(struct.pack('d', -1234567890123)) == gdal.CE_None
-        assert struct.unpack('d', var.Read()) == (-1234567890123, )
+        assert var.GetDataType().GetNumericDataType() == gdal.GDT_Int64
+        assert var.Write(struct.pack('q', -1234567890123)) == gdal.CE_None
+        assert struct.unpack('q', var.Read()) == (-1234567890123, )
 
         var = rg.CreateMDArray('var_as_nc_uint64', [],
                 gdal.ExtendedDataType.Create(gdal.GDT_Float64), ['NC_TYPE=NC_UINT64'])
-        assert var.GetDataType().GetNumericDataType() == gdal.GDT_Float64
-        assert var.Write(struct.pack('d', 1234567890123)) == gdal.CE_None
-        assert struct.unpack('d', var.Read()) == (1234567890123, )
+        assert var.GetDataType().GetNumericDataType() == gdal.GDT_UInt64
+        assert var.Write(struct.pack('Q', 1234567890123)) == gdal.CE_None
+        assert struct.unpack('Q', var.Read()) == (1234567890123, )
 
         # Test creation of compound data type
         comp0 = gdal.EDTComponent.Create('x', 0, gdal.ExtendedDataType.Create(gdal.GDT_Int16))
diff --git a/autotest/gdrivers/zarr_driver.py b/autotest/gdrivers/zarr_driver.py
index b1e164c76ff6..7de7ec02746a 100644
--- a/autotest/gdrivers/zarr_driver.py
+++ b/autotest/gdrivers/zarr_driver.py
@@ -52,15 +52,15 @@
                           [">i2", 'h', gdal.GDT_Int16, None, None],
                           ["<i4", 'i', gdal.GDT_Int32, None, None],
                           [">i4", 'i', gdal.GDT_Int32, None, None],
-                          ["<i8", 'q', gdal.GDT_Float64, None, None],
-                          [">i8", 'q', gdal.GDT_Float64, None, None],
+                          ["<i8", 'q', gdal.GDT_Int64, None, None],
+                          [">i8", 'q', gdal.GDT_Int64, None, None],
                           ["<u2", 'H', gdal.GDT_UInt16, None, None],
                           [">u2", 'H', gdal.GDT_UInt16, None, None],
                           ["<u4", 'I', gdal.GDT_UInt32, None, None],
                           [">u4", 'I', gdal.GDT_UInt32, None, None],
                           ["<u4", 'I', gdal.GDT_UInt32, 4000000000, 4000000000],
-                          ["<u8", 'Q', gdal.GDT_Float64, 4000000000, 4000000000],
-                          [">u8", 'Q', gdal.GDT_Float64, None, None],
+                          ["<u8", 'Q', gdal.GDT_UInt64, 4000000000, 4000000000],
+                          [">u8", 'Q', gdal.GDT_UInt64, None, None],
                           ["<f4", 'f', gdal.GDT_Float32, None, None],
                           [">f4", 'f', gdal.GDT_Float32, None, None],
                           ["<f4", 'f', gdal.GDT_Float32, 1.5, 1.5],
@@ -141,8 +141,6 @@ def test_zarr_basic(dtype, structtype, gdaltype, fill_value, nodata_value, use_o
 
         if structtype == 'b':
             structtype_read = 'h'
-        elif structtype in ('q', 'Q'):
-            structtype_read = 'd'
         else:
             structtype_read = structtype
 
diff --git a/autotest/utilities/test_gdalmdiminfo_lib.py b/autotest/utilities/test_gdalmdiminfo_lib.py
index 2ed0c8cbc791..c5deb909eecf 100755
--- a/autotest/utilities/test_gdalmdiminfo_lib.py
+++ b/autotest/utilities/test_gdalmdiminfo_lib.py
@@ -274,3 +274,33 @@ def test_gdalmdiminfo_lib_arrayoption():
     ret = gdal.MultiDimInfo('../gdrivers/data/netcdf/with_bounds.nc',
                             arrayoptions = ['SHOW_BOUNDS=NO'])
     assert len(ret['arrays']) == 1
+
+
+###############################################################################
+
+def test_gdalmdiminfo_lib_int64():
+
+    drv = gdal.GetDriverByName('MEM')
+    ds = drv.CreateMultiDimensional('')
+    rg = ds.GetRootGroup()
+    dim0 = rg.CreateDimension("dim0", "my_type", "my_direction", 1)
+    ar = rg.CreateMDArray("ar", [ dim0 ], gdal.ExtendedDataType.Create(gdal.GDT_Int64))
+    assert ar.Write(struct.pack('q', -10000000000)) == gdal.CE_None
+
+    ret = gdal.MultiDimInfo(ds, detailed = True)
+    assert ret['arrays']['ar']['values'] == [-10000000000]
+
+
+###############################################################################
+
+def test_gdalmdiminfo_lib_uint64():
+
+    drv = gdal.GetDriverByName('MEM')
+    ds = drv.CreateMultiDimensional('')
+    rg = ds.GetRootGroup()
+    dim0 = rg.CreateDimension("dim0", "my_type", "my_direction", 1)
+    ar = rg.CreateMDArray("ar", [ dim0 ], gdal.ExtendedDataType.Create(gdal.GDT_UInt64))
+    assert ar.Write(struct.pack('Q', 10000000000)) == gdal.CE_None
+
+    ret = gdal.MultiDimInfo(ds, detailed = True)
+    assert ret['arrays']['ar']['values'] == [10000000000]
diff --git a/data/gdalvrt.xsd b/data/gdalvrt.xsd
index 0f19f02dafc0..661e1d4d3d49 100644
--- a/data/gdalvrt.xsd
+++ b/data/gdalvrt.xsd
@@ -315,6 +315,8 @@
             <xs:enumeration value="Int16"/>
             <xs:enumeration value="UInt32"/>
             <xs:enumeration value="Int32"/>
+            <xs:enumeration value="UInt64"/>
+            <xs:enumeration value="Int64"/>
             <xs:enumeration value="Float32"/>
             <xs:enumeration value="Float64"/>
             <xs:enumeration value="CInt16"/>
diff --git a/doc/source/programs/gdal_create.rst b/doc/source/programs/gdal_create.rst
index 7fea7775f3cf..61570c7225ee 100644
--- a/doc/source/programs/gdal_create.rst
+++ b/doc/source/programs/gdal_create.rst
@@ -23,7 +23,7 @@ Synopsis
        [-outsize xsize ysize]
        [-bands count]
        [-burn value]*
-       [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
+       [-ot {Byte/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/
              CInt16/CInt32/CFloat32/CFloat64}] [-strict]
        [-a_srs srs_def] [-a_ullr ulx uly lrx lry] [-a_nodata value]
        [-mo "META-TAG=VALUE"]* [-q]
diff --git a/doc/source/programs/gdal_rasterize.rst b/doc/source/programs/gdal_rasterize.rst
index f870f5356fe0..3deac86b3fa4 100644
--- a/doc/source/programs/gdal_rasterize.rst
+++ b/doc/source/programs/gdal_rasterize.rst
@@ -21,7 +21,7 @@ Synopsis
         [-dialect dialect] [-of format] [-a_srs srs_def] [-to NAME=VALUE]*
         [-co "NAME=VALUE"]* [-a_nodata value] [-init value]*
         [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
-        [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
+        [-ot {Byte/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/
                 CInt16/CInt32/CFloat32/CFloat64}]
         [-optim {[AUTO]/VECTOR/RASTER}] [-q]
         <src_datasource> <dst_filename>
diff --git a/doc/source/programs/gdal_translate.rst b/doc/source/programs/gdal_translate.rst
index 0f4cda32b453..f833fc72d110 100644
--- a/doc/source/programs/gdal_translate.rst
+++ b/doc/source/programs/gdal_translate.rst
@@ -17,7 +17,7 @@ Synopsis
 
 
     gdal_translate [--help-general]
-        [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
+        [-ot {Byte/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/
                 CInt16/CInt32/CFloat32/CFloat64}] [-strict]
         [-if format]* [-of format]
         [-b band]* [-mask band] [-expand {gray|rgb|rgba}]
diff --git a/doc/source/user/raster_data_model.rst b/doc/source/user/raster_data_model.rst
index bb49df212321..c448d25e9093 100644
--- a/doc/source/user/raster_data_model.rst
+++ b/doc/source/user/raster_data_model.rst
@@ -196,7 +196,10 @@ A raster band is represented in GDAL with the :cpp:class:`GDALRasterBand` class.
 A raster band has the following properties:
 
 - A width and height in pixels and lines. This is the same as that defined for the dataset, if this is a full resolution band.
-- A datatype (GDALDataType). One of Byte, UInt16, Int16, UInt32, Int32, Float32, Float64, and the complex types CInt16, CInt32, CFloat32, and CFloat64.
+- A datatype (GDALDataType). One of Byte, UInt16, Int16, UInt32, Int32, UInt64, Int64, Float32, Float64, and the complex types CInt16, CInt32, CFloat32, and CFloat64.
+
+  UInt64 and Int64 data types have been added in GDAL 3.5. Beyond reading and write pixel values, their support is limited.  Some algorithms might use 64-bit floating-point internally (warping), as well as some methods like GetNoDataValue(), or even 32-bit floating point (overview, RasterIO resampling). Hence the range where exact values are preserved can be [0, 2^53] (or less if 32-bit floating-point is used).
+
 - A block size. This is a preferred (efficient) access chunk size. For tiled images this will be one tile. For scanline oriented images this will normally be one scanline.
 - A list of name/value pair metadata in the same format as the dataset, but of information that is potentially specific to this band.
 - An optional description string.
diff --git a/frmts/gtiff/cogdriver.cpp b/frmts/gtiff/cogdriver.cpp
index c40b8e04431e..8e565a185de5 100644
--- a/frmts/gtiff/cogdriver.cpp
+++ b/frmts/gtiff/cogdriver.cpp
@@ -1358,7 +1358,7 @@ void GDALRegister_COG()
     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/cog.html" );
 
     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
-                               "Byte UInt16 Int16 UInt32 Int32 Float32 "
+                               "Byte UInt16 Int16 UInt32 Int32 UInt64 Int64 Float32 "
                                "Float64 CInt16 CInt32 CFloat32 CFloat64" );
 
     poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
diff --git a/frmts/gtiff/geotiff.cpp b/frmts/gtiff/geotiff.cpp
index eb5f7b4021ca..29d2be091a5d 100644
--- a/frmts/gtiff/geotiff.cpp
+++ b/frmts/gtiff/geotiff.cpp
@@ -1431,6 +1431,10 @@ GTiffRasterBand::GTiffRasterBand( GTiffDataset *poDSIn, int nBandIn ):
             eDataType = GDT_CFloat32;
         else if( nSampleFormat == SAMPLEFORMAT_COMPLEXINT )
             eDataType = GDT_CInt32;
+        else if( nSampleFormat == SAMPLEFORMAT_INT )
+            eDataType = GDT_Int64;
+        else
+            eDataType = GDT_UInt64;
     }
     else if( nBitsPerSample == 128 )
     {
@@ -8527,6 +8531,18 @@ inline bool GTiffDataset::IsFirstPixelEqualToNoData( const void* pBuffer )
                *(static_cast<const GInt32*>(pBuffer)) ==
                         static_cast<GInt32>(dfEffectiveNoData);
     }
+    if( m_nBitsPerSample == 64 && eDT == GDT_UInt64 )
+    {
+        return GDALIsValueInRange<std::uint64_t>(dfEffectiveNoData) &&
+               *(static_cast<const std::uint64_t*>(pBuffer)) ==
+                        static_cast<std::uint64_t>(dfEffectiveNoData);
+    }
+    if( m_nBitsPerSample == 64 && eDT == GDT_Int64 )
+    {
+        return GDALIsValueInRange<std::int64_t>(dfEffectiveNoData) &&
+               *(static_cast<const std::int64_t*>(pBuffer)) ==
+                        static_cast<std::int64_t>(dfEffectiveNoData);
+    }
     if( m_nBitsPerSample == 32 && eDT == GDT_Float32 )
     {
         if( CPLIsNan(m_dfNoDataValue) )
@@ -9434,6 +9450,16 @@ template<> uint32_t AdjustValue<uint32_t>(uint32_t value, uint64_t nRoundUpBitTe
     return AdjustValueInt(value, nRoundUpBitTest);
 }
 
+template<> int64_t AdjustValue<int64_t>(int64_t value, uint64_t nRoundUpBitTest)
+{
+    return AdjustValueInt(value, nRoundUpBitTest);
+}
+
+template<> uint64_t AdjustValue<uint64_t>(uint64_t value, uint64_t nRoundUpBitTest)
+{
+    return AdjustValueInt(value, nRoundUpBitTest);
+}
+
 template<> float AdjustValue<float>(float value, uint64_t)
 {
     return std::nextafter(value, std::numeric_limits<float>::max());
@@ -9490,6 +9516,12 @@ template<> uint32_t RoundValueDiscardLsb<uint32_t, uint32_t>(const void* ptr,
     return RoundValueDiscardLsbUnsigned<uint32_t>(ptr, nMask, nRoundUpBitTest);
 }
 
+template<> uint64_t RoundValueDiscardLsb<uint64_t, uint64_t>(const void* ptr,
+                                                             uint64_t nMask,
+                                                             uint64_t nRoundUpBitTest)
+{
+    return RoundValueDiscardLsbUnsigned<uint64_t>(ptr, nMask, nRoundUpBitTest);
+}
 template<> int8_t RoundValueDiscardLsb<int8_t, int8_t>(const void* ptr,
                                                           uint64_t nMask,
                                                           uint64_t nRoundUpBitTest)
@@ -9511,6 +9543,13 @@ template<> int32_t RoundValueDiscardLsb<int32_t, int32_t>(const void* ptr,
     return RoundValueDiscardLsbSigned<int32_t>(ptr, nMask, nRoundUpBitTest);
 }
 
+template<> int64_t RoundValueDiscardLsb<int64_t, int64_t>(const void* ptr,
+                                                             uint64_t nMask,
+                                                             uint64_t nRoundUpBitTest)
+{
+    return RoundValueDiscardLsbSigned<int64_t>(ptr, nMask, nRoundUpBitTest);
+}
+
 template<> uint32_t RoundValueDiscardLsb<float, uint32_t>(const void* ptr,
                                                              uint64_t nMask,
                                                              uint64_t nRoundUpBitTest)
@@ -9784,6 +9823,44 @@ static void DiscardLsb( GByte* pabyBuffer, GPtrDiff_t nBytes, int iBand,
             pabyBuffer, nBytes, iBand, nBands, nPlanarConfig,
             panMaskOffsetLsb, bHasNoData, nNoDataValue);
     }
+    else if( nBitsPerSample == 64 && nSampleFormat == SAMPLEFORMAT_INT )
+    {
+        // FIXME: we should not rely on dfNoDataValue when we support native data type for nodata
+        int64_t nNoDataValue = 0;
+        if( bHasNoData &&
+            dfNoDataValue >= static_cast<double>(std::numeric_limits<int64_t>::min()) &&
+            dfNoDataValue <= static_cast<double>(std::numeric_limits<int64_t>::max()) &&
+            dfNoDataValue == static_cast<double>(static_cast<int64_t>(dfNoDataValue)) )
+        {
+            nNoDataValue = static_cast<int64_t>(dfNoDataValue);
+        }
+        else
+        {
+            bHasNoData = false;
+        }
+        DiscardLsbT<int64_t, int64_t>(
+            pabyBuffer, nBytes, iBand, nBands, nPlanarConfig,
+            panMaskOffsetLsb, bHasNoData, nNoDataValue);
+    }
+    else if( nBitsPerSample == 64 && nSampleFormat == SAMPLEFORMAT_UINT )
+    {
+        // FIXME: we should not rely on dfNoDataValue when we support native data type for nodata
+        uint64_t nNoDataValue = 0;
+        if( bHasNoData &&
+            dfNoDataValue >= static_cast<double>(std::numeric_limits<uint64_t>::min()) &&
+            dfNoDataValue <= static_cast<double>(std::numeric_limits<uint64_t>::max()) &&
+            dfNoDataValue == static_cast<double>(static_cast<uint64_t>(dfNoDataValue)) )
+        {
+            nNoDataValue = static_cast<uint64_t>(dfNoDataValue);
+        }
+        else
+        {
+            bHasNoData = false;
+        }
+        DiscardLsbT<uint64_t, uint64_t>(
+            pabyBuffer, nBytes, iBand, nBands, nPlanarConfig,
+            panMaskOffsetLsb, bHasNoData, nNoDataValue);
+    }
     else if( nBitsPerSample == 32 && nSampleFormat == SAMPLEFORMAT_IEEEFP )
     {
         float fNoDataValue = static_cast<float>(dfNoDataValue);
@@ -16545,7 +16622,7 @@ TIFF *GTiffDataset::CreateLL( const char * pszFilename,
 
     uint16_t l_nSampleFormat = 0;
     if( (eType == GDT_Byte && EQUAL(pszPixelType,"SIGNEDBYTE"))
-        || eType == GDT_Int16 || eType == GDT_Int32 )
+        || eType == GDT_Int16 || eType == GDT_Int32 || eType == GDT_Int64 )
         l_nSampleFormat = SAMPLEFORMAT_INT;
     else if( eType == GDT_CInt16 || eType == GDT_CInt32 )
         l_nSampleFormat = SAMPLEFORMAT_COMPLEXINT;
diff --git a/frmts/gtiff/gt_overview.cpp b/frmts/gtiff/gt_overview.cpp
index 1ee1f7a4bb42..86fc4d08aa5e 100644
--- a/frmts/gtiff/gt_overview.cpp
+++ b/frmts/gtiff/gt_overview.cpp
@@ -378,6 +378,16 @@ GTIFFBuildOverviewsEx( const char * pszFilename,
             nBandFormat = SAMPLEFORMAT_INT;
             break;
 
+          case GDT_UInt64:
+            nBandBits = 64;
+            nBandFormat = SAMPLEFORMAT_UINT;
+            break;
+
+          case GDT_Int64:
+            nBandBits = 64;
+            nBandFormat = SAMPLEFORMAT_INT;
+            break;
+
           case GDT_Float32:
             nBandBits = 32;
             nBandFormat = SAMPLEFORMAT_IEEEFP;
diff --git a/frmts/hdf4/hdf4dataset.cpp b/frmts/hdf4/hdf4dataset.cpp
index 8ec7c153c5db..e1f34e19920f 100644
--- a/frmts/hdf4/hdf4dataset.cpp
+++ b/frmts/hdf4/hdf4dataset.cpp
@@ -205,9 +205,9 @@ GDALDataType HDF4Dataset::GetDataType( int32 iNumType )
         case DFNT_UINT32:
             return GDT_UInt32;
         case DFNT_INT64:
-            return GDT_Unknown;
+            return GDT_Int64;
         case DFNT_UINT64:
-            return GDT_Unknown;
+            return GDT_UInt64;
         case DFNT_FLOAT32:
             return GDT_Float32;
         case DFNT_FLOAT64:
diff --git a/frmts/hdf4/hdf4imagedataset.cpp b/frmts/hdf4/hdf4imagedataset.cpp
index 11321b05191b..90034e3c4656 100644
--- a/frmts/hdf4/hdf4imagedataset.cpp
+++ b/frmts/hdf4/hdf4imagedataset.cpp
@@ -3949,6 +3949,38 @@ GDALDataset *HDF4ImageDataset::Create( const char * pszFilename,
     aiDimSizes[poDS->iYDim] = nYSize;
     aiDimSizes[poDS->iBandDim] = nBands;
 
+    const auto GetHDFType = [](GDALDataType eTypeIn)
+    {
+        switch ( eTypeIn )
+        {
+            case GDT_Float64:
+                return DFNT_FLOAT64;
+            case GDT_Float32:
+                return DFNT_FLOAT32;
+            case GDT_UInt64:
+                // SDCreate() doens't like it
+                return DFNT_UINT64;
+            case GDT_UInt32:
+                return DFNT_UINT32;
+            case GDT_UInt16:
+                return DFNT_UINT16;
+            case GDT_Int64:
+                // SDCreate() doens't like it
+                return DFNT_INT64;
+            case GDT_Int32:
+                return DFNT_INT32;
+            case GDT_Int16:
+                return DFNT_INT16;
+            case GDT_Byte:
+                return DFNT_UINT8;
+            default:
+                CPLError(CE_Warning, CPLE_NotSupported,
+                         "Datatype %s not supported. Defauting to Byte",
+                         GDALGetDataTypeName(eTypeIn));
+                return DFNT_UINT8;
+        }
+    };
+
     const char *pszSDSName = nullptr;
     int32 iSDS = -1;
 
@@ -3957,38 +3989,10 @@ GDALDataset *HDF4ImageDataset::Create( const char * pszFilename,
         for( int iBand = 0; iBand < nBands; iBand++ )
         {
             pszSDSName = CPLSPrintf( "Band%d", iBand );
-            switch ( eType )
-            {
-                case GDT_Float64:
-                    iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT64,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-                case GDT_Float32:
-                    iSDS = SDcreate( poDS-> hSD, pszSDSName, DFNT_FLOAT32,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-                case GDT_UInt32:
-                    iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT32,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-                case GDT_UInt16:
-                    iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT16,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-                case GDT_Int32:
-                    iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT32,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-                case GDT_Int16:
-                    iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT16,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-                case GDT_Byte:
-                default:
-                    iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT8,
-                                     poDS->iRank, aiDimSizes );
-                    break;
-            }
+            iSDS = SDcreate( poDS->hSD, pszSDSName, GetHDFType(eType),
+                             poDS->iRank, aiDimSizes );
+            if( iSDS < 0 )
+                break;
             SDendaccess( iSDS );
         }
     }
@@ -3996,38 +4000,8 @@ GDALDataset *HDF4ImageDataset::Create( const char * pszFilename,
     {
         pszSDSName = "3-dimensional Scientific Dataset";
         poDS->iDataset = 0;
-        switch ( eType )
-        {
-            case GDT_Float64:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT64,
-                                 poDS->iRank, aiDimSizes );
-                break;
-            case GDT_Float32:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT32,
-                                 poDS->iRank, aiDimSizes );
-                break;
-            case GDT_UInt32:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT32,
-                                 poDS->iRank, aiDimSizes );
-                break;
-            case GDT_UInt16:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT16,
-                                 poDS->iRank, aiDimSizes );
-                break;
-            case GDT_Int32:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT32,
-                                 poDS->iRank, aiDimSizes );
-                break;
-            case GDT_Int16:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT16,
-                                 poDS->iRank, aiDimSizes );
-                break;
-            case GDT_Byte:
-            default:
-                iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT8,
-                                 poDS->iRank, aiDimSizes );
-                break;
-        }
+        iSDS = SDcreate( poDS->hSD, pszSDSName, GetHDFType(eType),
+                         poDS->iRank, aiDimSizes );
     }
     else
     {
@@ -4089,6 +4063,7 @@ void GDALRegister_HDF4Image()
     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/hdf4.html" );
     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
                                "Byte Int16 UInt16 Int32 UInt32 "
+                               // "Int64 UInt64 "
                                "Float32 Float64" );
     poDriver->SetMetadataItem(
         GDAL_DMD_CREATIONOPTIONLIST,
diff --git a/frmts/hdf5/hdf5dataset.cpp b/frmts/hdf5/hdf5dataset.cpp
index 22afe724d82a..493758f1c212 100644
--- a/frmts/hdf5/hdf5dataset.cpp
+++ b/frmts/hdf5/hdf5dataset.cpp
@@ -180,6 +180,10 @@ GDALDataType HDF5Dataset::GetDataType(hid_t TypeID)
             return GDT_Int32;
         else if( H5Tequal(H5T_NATIVE_UINT,   TypeID) )
             return GDT_UInt32;
+        else if( H5Tequal(H5T_NATIVE_INT64,  TypeID) )
+            return GDT_Int64;
+        else if( H5Tequal(H5T_NATIVE_UINT64, TypeID) )
+            return GDT_UInt64;
         else if( H5Tequal(H5T_NATIVE_LONG,   TypeID) )
         {
 #if SIZEOF_UNSIGNED_LONG == 4
diff --git a/frmts/hdf5/hdf5multidim.cpp b/frmts/hdf5/hdf5multidim.cpp
index a6f67bcb8c69..2bbe6116b5b7 100644
--- a/frmts/hdf5/hdf5multidim.cpp
+++ b/frmts/hdf5/hdf5multidim.cpp
@@ -131,16 +131,6 @@ static GDALExtendedDataType BuildDataType(hid_t hDataType, bool& bHasVLen, bool&
         bNonNativeDataType = true;
         return GDALExtendedDataType::Create(GDT_Int16);
     }
-    else if( H5Tequal(H5T_NATIVE_LLONG, hDataType) )
-    {
-        bNonNativeDataType = true;
-        return GDALExtendedDataType::Create(GDT_Float64);
-    }
-    else if( H5Tequal(H5T_NATIVE_ULLONG, hDataType) )
-    {
-        bNonNativeDataType = true;
-        return GDALExtendedDataType::Create(GDT_Float64);
-    }
     else if( eDT != GDT_Unknown )
         return GDALExtendedDataType::Create(eDT);
     else if (klass == H5T_STRING )
@@ -1433,6 +1423,8 @@ static hid_t GetHDF5DataTypeFromGDALDataType(const GDALExtendedDataType& dt,
         case GDT_Int16: hBufferType = H5Tcopy(H5T_NATIVE_SHORT); break;
         case GDT_UInt32: hBufferType = H5Tcopy(H5T_NATIVE_UINT); break;
         case GDT_Int32: hBufferType = H5Tcopy(H5T_NATIVE_INT); break;
+        case GDT_UInt64: hBufferType = H5Tcopy(H5T_NATIVE_UINT64); break;
+        case GDT_Int64: hBufferType = H5Tcopy(H5T_NATIVE_INT64); break;
         case GDT_Float32: hBufferType = H5Tcopy(H5T_NATIVE_FLOAT); break;
         case GDT_Float64: hBufferType = H5Tcopy(H5T_NATIVE_DOUBLE); break;
         case GDT_CInt16:
@@ -1592,22 +1584,6 @@ static void CopyValue(const GByte* pabySrcBuffer, hid_t hSrcDataType,
                 &nVal, GDALExtendedDataType::Create(GDT_Int16),
                 pabyDstBuffer, dstDataType);
         }
-        else if( H5Tequal(H5T_NATIVE_LLONG, hSrcDataType) )
-        {
-            const double dfVal = static_cast<double>(
-                *reinterpret_cast<const GInt64*>(pabySrcBuffer));
-            GDALExtendedDataType::CopyValue(
-                &dfVal, GDALExtendedDataType::Create(GDT_Float64),
-                pabyDstBuffer, dstDataType);
-        }
-        else if( H5Tequal(H5T_NATIVE_ULLONG, hSrcDataType) )
-        {
-            const double dfVal = static_cast<double>(
-                *reinterpret_cast<const GUInt64*>(pabySrcBuffer));
-            GDALExtendedDataType::CopyValue(
-                &dfVal, GDALExtendedDataType::Create(GDT_Float64),
-                pabyDstBuffer, dstDataType);
-        }
         else
         {
             GDALDataType eDT = ::HDF5Dataset::GetDataType(hSrcDataType);
@@ -1744,7 +1720,9 @@ bool HDF5Array::IRead(const GUInt64* arrayStartIdx,
                 H5Tequal(hParent, H5T_NATIVE_USHORT) ||
                 H5Tequal(hParent, H5T_NATIVE_SHORT) ||
                 H5Tequal(hParent, H5T_NATIVE_UINT) ||
-                H5Tequal(hParent, H5T_NATIVE_INT) )
+                H5Tequal(hParent, H5T_NATIVE_INT) ||
+                H5Tequal(hParent, H5T_NATIVE_UINT64) ||
+                H5Tequal(hParent, H5T_NATIVE_INT64) )
             {
                 hBufferType = H5Tcopy(m_hNativeDT);
                 if( m_dt != bufferDataType )
@@ -2025,7 +2003,9 @@ bool HDF5Attribute::IRead(const GUInt64* arrayStartIdx,
                 H5Tequal(hParent, H5T_NATIVE_USHORT) ||
                 H5Tequal(hParent, H5T_NATIVE_SHORT) ||
                 H5Tequal(hParent, H5T_NATIVE_UINT) ||
-                H5Tequal(hParent, H5T_NATIVE_INT) )
+                H5Tequal(hParent, H5T_NATIVE_INT) ||
+                H5Tequal(hParent, H5T_NATIVE_UINT64) ||
+                H5Tequal(hParent, H5T_NATIVE_INT64))
             {
                 hBufferType = H5Tcopy(m_hNativeDT);
             }
diff --git a/frmts/mem/memdataset.cpp b/frmts/mem/memdataset.cpp
index c08c10eaccdc..48bc657e36d1 100644
--- a/frmts/mem/memdataset.cpp
+++ b/frmts/mem/memdataset.cpp
@@ -2471,7 +2471,7 @@ void GDALRegister_MEM()
     poDriver->SetMetadataItem( GDAL_DCAP_MULTIDIM_RASTER, "YES" );
     poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "In Memory Raster" );
     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
-                               "Byte Int16 UInt16 Int32 UInt32 Float32 Float64 "
+                               "Byte Int16 UInt16 Int32 UInt32 Int64 UInt64 Float32 Float64 "
                                "CInt16 CInt32 CFloat32 CFloat64" );
     poDriver->SetMetadataItem( GDAL_DCAP_COORDINATE_EPOCH, "YES" );
 
diff --git a/frmts/netcdf/netcdfdataset.cpp b/frmts/netcdf/netcdfdataset.cpp
index 340d72e9e988..27ee6451664e 100644
--- a/frmts/netcdf/netcdfdataset.cpp
+++ b/frmts/netcdf/netcdfdataset.cpp
@@ -391,6 +391,10 @@ netCDFRasterBand::netCDFRasterBand( const netCDFRasterBand::CONSTRUCTOR_OPEN&,
             eDataType = GDT_UInt16;
         else if( nc_datatype == NC_UINT )
             eDataType = GDT_UInt32;
+        else if( nc_datatype == NC_INT64 )
+            eDataType = GDT_Int64;
+        else if( nc_datatype == NC_UINT64 )
+            eDataType = GDT_UInt64;
 #endif
         else
         {
@@ -589,7 +593,8 @@ netCDFRasterBand::netCDFRasterBand( const netCDFRasterBand::CONSTRUCTOR_OPEN&,
 #ifdef NETCDF_HAS_NC4
     if( nc_datatype == NC_UBYTE ||
         nc_datatype == NC_USHORT ||
-        nc_datatype == NC_UINT )
+        nc_datatype == NC_UINT ||
+        nc_datatype == NC_UINT64 )
         bSignedData = false;
 #endif
 
@@ -779,6 +784,34 @@ netCDFRasterBand::netCDFRasterBand( const netCDFRasterBand::CONSTRUCTOR_CREATE&,
             nc_datatype = NC_DOUBLE;
             break;
 #ifdef NETCDF_HAS_NC4
+        case GDT_Int64:
+            if( poNCDFDS->eFormat == NCDF_FORMAT_NC4 )
+            {
+                nc_datatype = NC_INT64;
+            }
+            else
+            {
+                if( nBand == 1 )
+                    CPLError(CE_Warning, CPLE_AppDefined,
+                            "Unsupported GDAL datatype %s, treat as NC_DOUBLE.", "Int64");
+                nc_datatype = NC_DOUBLE;
+                eDataType = GDT_Float64;
+            }
+            break;
+        case GDT_UInt64:
+            if( poNCDFDS->eFormat == NCDF_FORMAT_NC4 )
+            {
+                nc_datatype = NC_UINT64;
+            }
+            else
+            {
+                if( nBand == 1 )
+                    CPLError(CE_Warning, CPLE_AppDefined,
+                            "Unsupported GDAL datatype %s, treat as NC_DOUBLE.", "UInt64");
+                nc_datatype = NC_DOUBLE;
+                eDataType = GDT_Float64;
+            }
+            break;
         case GDT_UInt16:
             if( poNCDFDS->eFormat == NCDF_FORMAT_NC4 )
             {
@@ -2004,6 +2037,22 @@ bool netCDFRasterBand::FetchNetcdfChunk( size_t xstart,
             CheckData<unsigned int>(pImage, pImageNC, edge[nBandXPos],
                                     nYChunkSize, false);
     }
+    else if( eDataType == GDT_Int64 )
+    {
+        status = nc_get_vara_longlong(cdfid, nZId, start, edge,
+                                      static_cast<long long *>(pImageNC));
+        if( status == NC_NOERR )
+            CheckData<std::int64_t>(pImage, pImageNC, edge[nBandXPos],
+                                    nYChunkSize, false);
+    }
+    else if( eDataType == GDT_UInt64 )
+    {
+        status = nc_get_vara_ulonglong(cdfid, nZId, start, edge,
+                                       static_cast<unsigned long long *>(pImageNC));
+        if( status == NC_NOERR )
+            CheckData<std::uint64_t>(pImage, pImageNC, edge[nBandXPos],
+                                     nYChunkSize, false);
+    }
     else if ( eDataType == GDT_CInt16 )
     {
         status = nc_get_vara(cdfid, nZId, start, edge,
@@ -2306,6 +2355,18 @@ CPLErr netCDFRasterBand::IWriteBlock( CPL_UNUSED int nBlockXOff,
         status = nc_put_vara_uint(cdfid, nZId, start, edge,
                                   static_cast<unsigned int *>(pImage));
     }
+    else if( eDataType == GDT_UInt64 &&
+             static_cast<netCDFDataset *>(poDS)->eFormat == NCDF_FORMAT_NC4 )
+    {
+        status = nc_put_vara_ulonglong(cdfid, nZId, start, edge,
+                                       static_cast<unsigned long long *>(pImage));
+    }
+    else if( eDataType == GDT_Int64 &&
+             static_cast<netCDFDataset *>(poDS)->eFormat == NCDF_FORMAT_NC4 )
+    {
+        status = nc_put_vara_longlong(cdfid, nZId, start, edge,
+                                      static_cast<long long *>(pImage));
+    }
 #endif
     else
     {
@@ -8948,9 +9009,19 @@ netCDFDataset::Create( const char *pszFilename,
 
     CPLMutexHolderD(&hNCMutex);
 
+    CPLStringList aosOptions(CSLDuplicate(papszOptions));
+#ifdef NETCDF_HAS_NC4
+    if( aosOptions.FetchNameValue("FORMAT") == nullptr &&
+        (eType == GDT_UInt16 || eType == GDT_UInt32 ||
+         eType == GDT_UInt64 || eType == GDT_Int64) )
+    {
+        CPLDebug("netCDF", "Selecting FORMAT=NC4 due to data type");
+        aosOptions.SetNameValue("FORMAT", "NC4");
+    }
+#endif
     netCDFDataset *poDS = netCDFDataset::CreateLL(pszFilename,
                                                   nXSize, nYSize, nBands,
-                                                  papszOptions);
+                                                  aosOptions.List());
 
     if( !poDS )
         return nullptr;
@@ -9088,7 +9159,7 @@ netCDFDataset::CreateCopy( const char *pszFilename, GDALDataset *poSrcDS,
         return nullptr;
     }
 
-    GDALDataType eDT;
+    GDALDataType eDT = GDT_Unknown;
     GDALRasterBand *poSrcBand = nullptr;
     for( int iBand = 1; iBand <= nBands; iBand++ )
     {
@@ -9107,9 +9178,19 @@ netCDFDataset::CreateCopy( const char *pszFilename, GDALDataset *poSrcDS,
         return nullptr;
 
     // Same as in Create().
+    CPLStringList aosOptions(CSLDuplicate(papszOptions));
+#ifdef NETCDF_HAS_NC4
+    if( aosOptions.FetchNameValue("FORMAT") == nullptr &&
+        (eDT == GDT_UInt16 || eDT == GDT_UInt32 ||
+         eDT == GDT_UInt64 || eDT == GDT_Int64) )
+    {
+        CPLDebug("netCDF", "Selecting FORMAT=NC4 due to data type");
+        aosOptions.SetNameValue("FORMAT", "NC4");
+    }
+#endif
     netCDFDataset *poDS = netCDFDataset::CreateLL(pszFilename,
                                                    nXSize, nYSize, nBands,
-                                                   papszOptions);
+                                                   aosOptions.List());
     if( !poDS )
         return nullptr;
 
@@ -9399,17 +9480,41 @@ netCDFDataset::CreateCopy( const char *pszFilename, GDALDataset *poSrcDS,
             eErr = NCDFCopyBand<GByte>(poSrcBand, poDstBand, nXSize, nYSize,
                                        GDALScaledProgress, pScaledProgress);
         }
-        else if( (eDT == GDT_UInt16) || (eDT == GDT_Int16) )
+        else if( eDT == GDT_UInt16 )
         {
-            CPLDebug("GDAL_netCDF", "GInt16 Band#%d", iBand);
+            CPLDebug("GDAL_netCDF", "GUInt16 Band#%d", iBand);
             eErr = NCDFCopyBand<GInt16>(poSrcBand, poDstBand, nXSize, nYSize,
                                         GDALScaledProgress, pScaledProgress);
         }
-        else if( (eDT == GDT_UInt32) || (eDT == GDT_Int32) )
+        else if( eDT == GDT_Int16 )
         {
             CPLDebug("GDAL_netCDF", "GInt16 Band#%d", iBand);
+            eErr = NCDFCopyBand<GUInt16>(poSrcBand, poDstBand, nXSize, nYSize,
+                                         GDALScaledProgress, pScaledProgress);
+        }
+        else if( eDT == GDT_UInt32 )
+        {
+            CPLDebug("GDAL_netCDF", "GUInt32 Band#%d", iBand);
+            eErr = NCDFCopyBand<GUInt32>(poSrcBand, poDstBand, nXSize, nYSize,
+                                         GDALScaledProgress, pScaledProgress);
+        }
+        else if( eDT == GDT_Int32 )
+        {
+            CPLDebug("GDAL_netCDF", "GInt32 Band#%d", iBand);
             eErr = NCDFCopyBand<GInt32>(poSrcBand, poDstBand, nXSize, nYSize,
-                                        GDALScaledProgress, pScaledProgress);
+                                         GDALScaledProgress, pScaledProgress);
+        }
+        else if( eDT == GDT_UInt64 )
+        {
+            CPLDebug("GDAL_netCDF", "GUInt64 Band#%d", iBand);
+            eErr = NCDFCopyBand<std::uint64_t>(poSrcBand, poDstBand, nXSize, nYSize,
+                                               GDALScaledProgress, pScaledProgress);
+        }
+        else if( eDT == GDT_Int64 )
+        {
+            CPLDebug("GDAL_netCDF", "GInt64 Band#%d", iBand);
+            eErr = NCDFCopyBand<std::int64_t>(poSrcBand, poDstBand, nXSize, nYSize,
+                                              GDALScaledProgress, pScaledProgress);
         }
         else if( eDT == GDT_Float32 )
         {
@@ -9728,6 +9833,14 @@ void GDALRegister_netCDF()
     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Network Common Data Format");
     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/netcdf.html");
     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "nc");
+    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
+#ifdef NETCDF_HAS_NC4
+                              "Byte UInt16 Int16 UInt32 Int32 Int64 UInt64 "
+#else
+                              "Byte Int16 Int32 "
+#endif
+                              "Float32 Float64 "
+                              "CInt16 CInt32 CFloat32 CFloat64" );
     poDriver->SetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST,
 "<CreationOptionList>"
 "   <Option name='FORMAT' type='string-select' default='NC'>"
diff --git a/frmts/netcdf/netcdfmultidim.cpp b/frmts/netcdf/netcdfmultidim.cpp
index 3a1dac8cad95..bf209088e0fd 100644
--- a/frmts/netcdf/netcdfmultidim.cpp
+++ b/frmts/netcdf/netcdfmultidim.cpp
@@ -708,6 +708,8 @@ static int CreateOrGetType(int gid, const GDALExtendedDataType& oType)
             case GDT_Int16: nTypeId = NC_SHORT; break;
             case GDT_UInt32: nTypeId = NC_UINT; break;
             case GDT_Int32: nTypeId = NC_INT; break;
+            case GDT_UInt64: nTypeId = NC_UINT64; break;
+            case GDT_Int64: nTypeId = NC_INT64; break;
             case GDT_Float32: nTypeId = NC_FLOAT; break;
             case GDT_Float64: nTypeId = NC_DOUBLE; break;
             case GDT_CInt16:
@@ -1915,9 +1917,15 @@ static bool BuildDataType(int gid, int varid, int nVarType,
             eDataType = GDT_UInt32;
         }
         else if( nVarType == NC_INT64 )
-            eDataType = GDT_Float64; // approximation
+        {
+            bPerfectDataTypeMatch = true;
+            eDataType = GDT_Int64;
+        }
         else if( nVarType == NC_UINT64 )
-            eDataType = GDT_Float64; // approximation
+        {
+            bPerfectDataTypeMatch = true;
+            eDataType = GDT_UInt64;
+        }
         else
         {
             CPLError(CE_Failure, CPLE_NotSupported,
diff --git a/frmts/vrt/pixelfunctions.cpp b/frmts/vrt/pixelfunctions.cpp
index 940bbe2fca88..e24aa10ba70c 100644
--- a/frmts/vrt/pixelfunctions.cpp
+++ b/frmts/vrt/pixelfunctions.cpp
@@ -45,6 +45,9 @@ template<typename T> inline double GetSrcVal(const void* pSource, GDALDataType e
         case GDT_Int16: return static_cast<const GInt16*>(pSource)[ii];
         case GDT_UInt32: return static_cast<const GUInt32*>(pSource)[ii];
         case GDT_Int32: return static_cast<const GInt32*>(pSource)[ii];
+        // Precision loss currently for int64/uint64
+        case GDT_UInt64: return static_cast<double>(static_cast<const uint64_t*>(pSource)[ii]);
+        case GDT_Int64: return static_cast<double>(static_cast<const int64_t*>(pSource)[ii]);
         case GDT_Float32: return static_cast<const float*>(pSource)[ii];
         case GDT_Float64: return static_cast<const double*>(pSource)[ii];
         case GDT_CInt16: return static_cast<const GInt16*>(pSource)[2 * ii];
diff --git a/frmts/vrt/vrtdriver.cpp b/frmts/vrt/vrtdriver.cpp
index a69f530a4d1d..7383b86a943d 100644
--- a/frmts/vrt/vrtdriver.cpp
+++ b/frmts/vrt/vrtdriver.cpp
@@ -439,7 +439,8 @@ void GDALRegister_VRT()
     poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "vrt" );
     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/vrt.html" );
     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
-                               "Byte Int16 UInt16 Int32 UInt32 Float32 Float64 "
+                               "Byte Int16 UInt16 Int32 UInt32 Int64 UInt64 "
+                               "Float32 Float64 "
                                "CInt16 CInt32 CFloat32 CFloat64" );
 
     poDriver->pfnOpen = VRTDataset::Open;
diff --git a/frmts/zarr/zarr_array.cpp b/frmts/zarr/zarr_array.cpp
index deb7ed789377..7bf6a263b429 100644
--- a/frmts/zarr/zarr_array.cpp
+++ b/frmts/zarr/zarr_array.cpp
@@ -456,20 +456,6 @@ static void EncodeElt(const std::vector<DtypeElt>& elts,
                         CPL_SWAP32(*reinterpret_cast<const uint32_t*>(pSrc + elt.gdalOffset + 4));
                     memcpy(pDst + elt.nativeOffset + 4, &val, sizeof(val));
                 }
-                else if( elt.nativeType == DtypeElt::NativeType::SIGNED_INT )
-                {
-                    CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                    const double dbl = *reinterpret_cast<const double*>(pSrc + elt.gdalOffset);
-                    int64_t val = CPL_SWAP64(static_cast<int64_t>(dbl));
-                    memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
-                }
-                else if( elt.nativeType == DtypeElt::NativeType::UNSIGNED_INT )
-                {
-                    CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                    const double dbl = *reinterpret_cast<const double*>(pSrc + elt.gdalOffset);
-                    uint64_t val = CPL_SWAP64(static_cast<uint64_t>(dbl));
-                    memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
-                }
                 else
                 {
                     const uint64_t val = CPL_SWAP64(
@@ -511,22 +497,6 @@ static void EncodeElt(const std::vector<DtypeElt>& elts,
                 const uint16_t uint16Val = CPLFloatToHalf(uint32Val, bHasWarned);
                 memcpy(pDst + elt.nativeOffset, &uint16Val, sizeof(uint16Val));
             }
-            else if( elt.nativeType == DtypeElt::NativeType::SIGNED_INT &&
-                     elt.nativeSize == 8 )
-            {
-                CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                const double dbl = *reinterpret_cast<const double*>(pSrc + elt.gdalOffset);
-                const int64_t val = static_cast<int64_t>(dbl);
-                memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
-            }
-            else if( elt.nativeType == DtypeElt::NativeType::UNSIGNED_INT &&
-                     elt.nativeSize == 8 )
-            {
-                CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                const double dbl = *reinterpret_cast<const double*>(pSrc + elt.gdalOffset);
-                const uint64_t val = static_cast<uint64_t>(dbl);
-                memcpy(pDst + elt.nativeOffset, &val, sizeof(val));
-            }
             else
             {
                 CPLAssert(false);
@@ -1101,22 +1071,6 @@ static void DecodeSourceElt(const std::vector<DtypeElt>& elts,
                     *reinterpret_cast<uint32_t*>(pDst + elt.gdalOffset + 4) =
                         CPL_SWAP32(val);
                 }
-                else if( elt.nativeType == DtypeElt::NativeType::SIGNED_INT )
-                {
-                    CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                    int64_t val;
-                    memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
-                    *reinterpret_cast<double*>(pDst + elt.gdalOffset) = static_cast<double>(
-                        CPL_SWAP64(val));
-                }
-                else if( elt.nativeType == DtypeElt::NativeType::UNSIGNED_INT )
-                {
-                    CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                    uint64_t val;
-                    memcpy(&val, pSrc + elt.nativeOffset, sizeof(val));
-                    *reinterpret_cast<double*>(pDst + elt.gdalOffset) = static_cast<double>(
-                        CPL_SWAP64(val));
-                }
                 else
                 {
                     uint64_t val;
@@ -1158,24 +1112,6 @@ static void DecodeSourceElt(const std::vector<DtypeElt>& elts,
                 uint32_t uint32Val = CPLHalfToFloat(uint16Val);
                 memcpy(pDst + elt.gdalOffset, &uint32Val, sizeof(uint32Val));
             }
-            else if( elt.nativeType == DtypeElt::NativeType::SIGNED_INT &&
-                     elt.nativeSize == 8 )
-            {
-                CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                int64_t intVal;
-                memcpy(&intVal, pSrc + elt.nativeOffset, sizeof(intVal));
-                double dblVal = static_cast<double>(intVal);
-                memcpy(pDst + elt.gdalOffset, &dblVal, sizeof(dblVal));
-            }
-            else if( elt.nativeType == DtypeElt::NativeType::UNSIGNED_INT &&
-                     elt.nativeSize == 8 )
-            {
-                CPLAssert( elt.gdalType.GetNumericDataType() == GDT_Float64 );
-                uint64_t intVal;
-                memcpy(&intVal, pSrc + elt.nativeOffset, sizeof(intVal));
-                double dblVal = static_cast<double>(intVal);
-                memcpy(pDst + elt.gdalOffset, &dblVal, sizeof(dblVal));
-            }
             else
             {
                 CPLAssert(false);
@@ -2988,8 +2924,7 @@ static GDALExtendedDataType ParseDtype(bool isZarrV2,
             else if( chType == 'i' && nBytes == 8 )
             {
                 elt.nativeType = DtypeElt::NativeType::SIGNED_INT;
-                elt.gdalTypeIsApproxOfNative = true;
-                eDT = GDT_Float64;
+                eDT = GDT_Int64;
             }
             else if( chType == 'u' && nBytes == 2 )
             {
@@ -3004,8 +2939,7 @@ static GDALExtendedDataType ParseDtype(bool isZarrV2,
             else if( chType == 'u' && nBytes == 8 )
             {
                 elt.nativeType = DtypeElt::NativeType::UNSIGNED_INT;
-                elt.gdalTypeIsApproxOfNative = true;
-                eDT = GDT_Float64;
+                eDT = GDT_UInt64;
             }
             else if( chType == 'f' && nBytes == 2 )
             {
diff --git a/frmts/zarr/zarr_group.cpp b/frmts/zarr/zarr_group.cpp
index 977b4b6a65b0..3cbfa88343db 100644
--- a/frmts/zarr/zarr_group.cpp
+++ b/frmts/zarr/zarr_group.cpp
@@ -912,6 +912,18 @@ static CPLJSONObject FillDTypeElts(const GDALExtendedDataType& oDataType,
                     dtype.Set(dummy, "<i4");
                     break;
                 }
+                case GDT_UInt64:
+                {
+                    elt.nativeType = DtypeElt::NativeType::UNSIGNED_INT;
+                    dtype.Set(dummy, "<u8");
+                    break;
+                }
+                case GDT_Int64:
+                {
+                    elt.nativeType = DtypeElt::NativeType::SIGNED_INT;
+                    dtype.Set(dummy, "<i8");
+                    break;
+                }
                 case GDT_Float32:
                 {
                     elt.nativeType = DtypeElt::NativeType::IEEEFP;
@@ -955,7 +967,7 @@ static CPLJSONObject FillDTypeElts(const GDALExtendedDataType& oDataType,
                 }
                 case GDT_TypeCount:
                 {
-                    static_assert(GDT_TypeCount == GDT_CFloat64 + 1, "GDT_TypeCount == GDT_CFloat64 + 1");
+                    static_assert(GDT_TypeCount == GDT_Int64 + 1, "GDT_TypeCount == GDT_Int64 + 1");
                     break;
                 }
             }
@@ -1242,6 +1254,8 @@ std::shared_ptr<GDALMDArray> ZarrGroupV2::CreateMDArray(
                 case GDT_Int16: oFilter.Add("dtype", "<i2"); break;
                 case GDT_UInt32: oFilter.Add("dtype", "<u4"); break;
                 case GDT_Int32: oFilter.Add("dtype", "<i4"); break;
+                case GDT_UInt64: oFilter.Add("dtype", "<u8"); break;
+                case GDT_Int64: oFilter.Add("dtype", "<i8"); break;
                 case GDT_Float32: oFilter.Add("dtype", "<f4"); break;
                 case GDT_Float64: oFilter.Add("dtype", "<f8"); break;
                 case GDT_CInt16: oFilter.Add("dtype", "<i2"); break;
diff --git a/frmts/zarr/zarrdriver.cpp b/frmts/zarr/zarrdriver.cpp
index 1b9f2b33a512..72a4264b6474 100644
--- a/frmts/zarr/zarrdriver.cpp
+++ b/frmts/zarr/zarrdriver.cpp
@@ -1259,7 +1259,7 @@ void GDALRegister_Zarr()
     poDriver->SetMetadataItem( GDAL_DCAP_MULTIDIM_RASTER, "YES" );
     poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Zarr" );
     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
-                               "Byte Int16 UInt16 Int32 UInt32 "
+                               "Byte Int16 UInt16 Int32 UInt32 Int64 UInt64 "
                                "Float32 Float64 CFloat32 CFloat64" );
     poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
     poDriver->SetMetadataItem( GDAL_DMD_SUBDATASETS, "YES" );
diff --git a/gcore/gdal.h b/gcore/gdal.h
index 58f566dc4df9..6c5237bb6e03 100644
--- a/gcore/gdal.h
+++ b/gcore/gdal.h
@@ -66,8 +66,8 @@ typedef enum {
     /*! Sixteen bit signed integer */           GDT_Int16 = 3,
     /*! Thirty two bit unsigned integer */      GDT_UInt32 = 4,
     /*! Thirty two bit signed integer */        GDT_Int32 = 5,
-    /* TODO?(#6879): GDT_UInt64 */
-    /* TODO?(#6879): GDT_Int64 */
+    /*! 64 bit unsigned integer (GDAL >= 3.5)*/ GDT_UInt64 = 12,
+    /*! 64 bit signed integer  (GDAL >= 3.5)*/  GDT_Int64 = 13,
     /*! Thirty two bit floating point */        GDT_Float32 = 6,
     /*! Sixty four bit floating point */        GDT_Float64 = 7,
     /*! Complex Int16 */                        GDT_CInt16 = 8,
@@ -75,7 +75,7 @@ typedef enum {
     /* TODO?(#6879): GDT_CInt64 */
     /*! Complex Float32 */                      GDT_CFloat32 = 10,
     /*! Complex Float64 */                      GDT_CFloat64 = 11,
-    GDT_TypeCount = 12          /* maximum type # + 1 */
+    GDT_TypeCount = 14          /* maximum type # + 1 */
 } GDALDataType;
 
 int CPL_DLL CPL_STDCALL GDALGetDataTypeSize( GDALDataType );  // Deprecated.
diff --git a/gcore/gdal_misc.cpp b/gcore/gdal_misc.cpp
index 6211be602dcd..1a9f95c8ba93 100644
--- a/gcore/gdal_misc.cpp
+++ b/gcore/gdal_misc.cpp
@@ -111,6 +111,8 @@ static int GetDataTypeElementSizeBits( GDALDataType eDataType )
 
       case GDT_Float64:
       case GDT_CFloat64:
+      case GDT_UInt64:
+      case GDT_Int64:
         return 64;
 
       default:
@@ -214,6 +216,10 @@ static int GetMinBitsForValue(double dValue)
         if( dValue <= std::numeric_limits<GUInt32>::max() &&
             dValue >= std::numeric_limits<GUInt32>::min() )
             return 32;
+
+        if( dValue <= static_cast<double>(std::numeric_limits<std::uint64_t>::max()) &&
+            dValue >= static_cast<double>(std::numeric_limits<std::uint64_t>::min()) )
+            return 64;
     }
     else if( static_cast<float>(dValue) == dValue )
     {
@@ -268,6 +274,9 @@ GDALDataType CPL_STDCALL GDALFindDataType(
         return GDT_UInt32;
     }
 
+    if( nBits == 64 && !bFloating && !bComplex )
+        return bSigned ? GDT_Int64 : GDT_UInt64;
+
     if( bComplex )
         return GDT_CFloat64;
 
@@ -332,6 +341,8 @@ int CPL_STDCALL GDALGetDataTypeSizeBytes( GDALDataType eDataType )
       case GDT_Float64:
       case GDT_CInt32:
       case GDT_CFloat32:
+      case GDT_UInt64:
+      case GDT_Int64:
         return 8;
 
       case GDT_CFloat64:
@@ -463,6 +474,8 @@ int CPL_STDCALL GDALDataTypeIsInteger( GDALDataType eDataType )
       case GDT_UInt32:
       case GDT_CInt16:
       case GDT_CInt32:
+      case GDT_UInt64:
+      case GDT_Int64:
         return TRUE;
 
       default:
@@ -488,6 +501,7 @@ int CPL_STDCALL GDALDataTypeIsSigned( GDALDataType eDataType )
       case GDT_Byte:
       case GDT_UInt16:
       case GDT_UInt32:
+      case GDT_UInt64:
         return FALSE;
 
       default:
@@ -545,11 +559,19 @@ int CPL_STDCALL GDALDataTypeIsConversionLossy( GDALDataType eTypeFrom,
 
     if( eTypeTo == GDT_Float32 && (eTypeFrom == GDT_Int32 ||
                                    eTypeFrom == GDT_UInt32 ||
+                                   eTypeFrom == GDT_Int64 ||
+                                   eTypeFrom == GDT_UInt64 ||
                                    eTypeFrom == GDT_Float64) )
     {
         return TRUE;
     }
 
+    if( eTypeTo == GDT_Float64 && (eTypeFrom == GDT_Int64 ||
+                                   eTypeFrom == GDT_UInt64) )
+    {
+        return TRUE;
+    }
+
     return FALSE;
 }
 
@@ -594,6 +616,12 @@ const char * CPL_STDCALL GDALGetDataTypeName( GDALDataType eDataType )
       case GDT_Int32:
         return "Int32";
 
+      case GDT_UInt64:
+        return "UInt64";
+
+      case GDT_Int64:
+        return "Int64";
+
       case GDT_Float32:
         return "Float32";
 
@@ -659,12 +687,12 @@ template<class T> static inline void ClampAndRound(
 {
     // TODO(schwehr): Rework this template.  ::min() versus ::lowest.
 
-    if (dfValue < std::numeric_limits<T>::min())
+    if (dfValue < static_cast<double>(std::numeric_limits<T>::min()))
     {
         bClamped = true;
         dfValue = static_cast<double>(std::numeric_limits<T>::min());
     }
-    else if (dfValue > std::numeric_limits<T>::max())
+    else if (dfValue > static_cast<double>(std::numeric_limits<T>::max()))
     {
         bClamped = true;
         dfValue = static_cast<double>(std::numeric_limits<T>::max());
@@ -713,6 +741,12 @@ double GDALAdjustValueToDataType(
         case GDT_UInt32:
             ClampAndRound<GUInt32>(dfValue, bClamped, bRounded);
             break;
+        case GDT_Int64:
+            ClampAndRound<std::int64_t>(dfValue, bClamped, bRounded);
+            break;
+        case GDT_UInt64:
+            ClampAndRound<std::uint64_t>(dfValue, bClamped, bRounded);
+            break;
         case GDT_Float32:
         {
             if( !CPLIsFinite(dfValue) )
@@ -1128,6 +1162,12 @@ GDALGetRandomRasterSample( GDALRasterBandH hBand, int nSamples,
                   case GDT_Int32:
                     dfValue = reinterpret_cast<const GInt32*>(pDataRef)[iOffset];
                     break;
+                  case GDT_UInt64:
+                    dfValue = static_cast<double>(reinterpret_cast<const std::uint64_t*>(pDataRef)[iOffset]);
+                    break;
+                  case GDT_Int64:
+                    dfValue = static_cast<double>(reinterpret_cast<const std::int64_t*>(pDataRef)[iOffset]);
+                    break;
                   case GDT_Float32:
                     dfValue = reinterpret_cast<const float*>(pDataRef)[iOffset];
                     break;
diff --git a/gcore/gdal_priv_templates.hpp b/gcore/gdal_priv_templates.hpp
index aeadec1b7842..7ae2a1d3cde4 100644
--- a/gcore/gdal_priv_templates.hpp
+++ b/gcore/gdal_priv_templates.hpp
@@ -32,6 +32,7 @@
 
 #include "cpl_port.h"
 
+#include <cstdint>
 #include <limits>
 
 /************************************************************************/
@@ -112,8 +113,8 @@ inline T GDALClampValue(const T tValue, const T tMax, const T tMin)
  */
 template <class T> inline bool GDALIsValueInRange(double dfValue)
 {
-    return dfValue >= std::numeric_limits<T>::min() &&
-           dfValue <= std::numeric_limits<T>::max();
+    return dfValue >= static_cast<double>(std::numeric_limits<T>::min()) &&
+           dfValue <= static_cast<double>(std::numeric_limits<T>::max());
 }
 
 template <> inline bool GDALIsValueInRange<double>(double dfValue)
@@ -154,7 +155,7 @@ template<class Tin> struct sGDALCopyWord<Tin, double>
 {
     static inline void f(const Tin tValueIn, double &dfValueOut)
     {
-        dfValueOut = tValueIn;
+        dfValueOut = static_cast<double>(tValueIn);
     }
 };
 
@@ -269,6 +270,24 @@ template<> struct sGDALCopyWord<double, int>
     }
 };
 
+template<> struct sGDALCopyWord<double, std::int64_t>
+{
+    static inline void f(const double dfValueIn, std::int64_t &nValueOut)
+    {
+        if( CPLIsNan(dfValueIn) )
+        {
+            nValueOut = 0;
+            return;
+        }
+        double dfMaxVal, dfMinVal;
+        GDALGetDataLimits<double, std::int64_t>(dfMaxVal, dfMinVal);
+        double dfValue = dfValueIn >= 0.0 ? dfValueIn + 0.5 :
+            dfValueIn - 0.5;
+        nValueOut = static_cast<std::int64_t>(
+            GDALClampValue(dfValue, dfMaxVal, dfMinVal));
+    }
+};
+
 template<> struct sGDALCopyWord<double, short>
 {
     static inline void f(const double dfValueIn, short &nValueOut)
@@ -330,6 +349,49 @@ template<> struct sGDALCopyWord<float, unsigned int>
     }
 };
 
+// Roundoff occurs for Float32 -> std::int64_t for max/min. Overload GDALCopyWord
+// specifically for this case.
+template<> struct sGDALCopyWord<float, std::int64_t>
+{
+    static inline void f(const float fValueIn, std::int64_t &nValueOut)
+    {
+        if (fValueIn >= static_cast<float>(std::numeric_limits<std::int64_t>::max()))
+        {
+            nValueOut = std::numeric_limits<std::int64_t>::max();
+        }
+        else if (fValueIn <= static_cast<float>(std::numeric_limits<std::int64_t>::min()))
+        {
+            nValueOut = std::numeric_limits<std::int64_t>::min();
+        }
+        else
+        {
+            nValueOut = static_cast<std::int64_t>(fValueIn > 0.0f ?
+                fValueIn + 0.5f : fValueIn - 0.5f);
+        }
+    }
+};
+
+// Roundoff occurs for Float32 -> std::uint64_t for max. Overload GDALCopyWord
+// specifically for this case.
+template<> struct sGDALCopyWord<float, std::uint64_t>
+{
+    static inline void f(const float fValueIn, std::uint64_t &nValueOut)
+    {
+        if (fValueIn >= static_cast<float>(std::numeric_limits<std::uint64_t>::max()))
+        {
+            nValueOut = std::numeric_limits<std::uint64_t>::max();
+        }
+        else if (fValueIn <= static_cast<float>(std::numeric_limits<std::uint64_t>::min()))
+        {
+            nValueOut = std::numeric_limits<std::uint64_t>::min();
+        }
+        else
+        {
+            nValueOut = static_cast<std::uint64_t>(fValueIn + 0.5f);
+        }
+    }
+};
+
 /**
  * Copy a single word, optionally rounding if appropriate (i.e. going
  * from the float to the integer case). Note that this is the function
diff --git a/gcore/gdalmultidim.cpp b/gcore/gdalmultidim.cpp
index 9314160be1cf..5601cd601596 100644
--- a/gcore/gdalmultidim.cpp
+++ b/gcore/gdalmultidim.cpp
@@ -819,7 +819,9 @@ bool GDALGroup::CopyFrom( const std::shared_ptr<GDALGroup>& poDstRootGroup,
                                 eAutoScaleType != GDT_UInt16 &&
                                 eAutoScaleType != GDT_Int16 &&
                                 eAutoScaleType != GDT_UInt32 &&
-                                eAutoScaleType != GDT_Int32 )
+                                eAutoScaleType != GDT_Int32 &&
+                                eAutoScaleType != GDT_UInt64 &&
+                                eAutoScaleType != GDT_Int64 )
                             {
                                 CPLError(CE_Failure, CPLE_NotSupported,
                                          "Unsupported value for AUTOSCALE_DATA_TYPE");
@@ -869,7 +871,8 @@ bool GDALGroup::CopyFrom( const std::shared_ptr<GDALGroup>& poDstRootGroup,
                 double dfDTMin = 0;
                 double dfDTMax = 0;
 #define setDTMinMax(ctype) do \
-    { dfDTMin = std::numeric_limits<ctype>::min(); dfDTMax = std::numeric_limits<ctype>::max(); } while(0)
+    { dfDTMin = static_cast<double>(std::numeric_limits<ctype>::min()); \
+      dfDTMax = static_cast<double>(std::numeric_limits<ctype>::max()); } while(0)
 
                 switch( eAutoScaleType )
                 {
@@ -878,6 +881,8 @@ bool GDALGroup::CopyFrom( const std::shared_ptr<GDALGroup>& poDstRootGroup,
                     case GDT_Int16:  setDTMinMax(GInt16); break;
                     case GDT_UInt32: setDTMinMax(GUInt32); break;
                     case GDT_Int32:  setDTMinMax(GInt32); break;
+                    case GDT_UInt64: setDTMinMax(std::uint64_t); break;
+                    case GDT_Int64:  setDTMinMax(std::int64_t); break;
                     default:
                         CPLAssert(false);
                 }
@@ -1340,6 +1345,12 @@ bool GDALExtendedDataType::CopyValue(const void* pSrc,
             case GDT_Int32:
                 str = CPLSPrintf("%d", *static_cast<const GInt32*>(pSrc));
                 break;
+            case GDT_UInt64:
+                str = CPLSPrintf(CPL_FRMT_GUIB, static_cast<GUIntBig>(*static_cast<const std::uint64_t*>(pSrc)));
+                break;
+            case GDT_Int64:
+                str = CPLSPrintf(CPL_FRMT_GIB, static_cast<GIntBig>(*static_cast<const std::int64_t*>(pSrc)));
+                break;
             case GDT_Float32:
                 str = CPLSPrintf("%.9g", *static_cast<const float*>(pSrc));
                 break;
@@ -1383,10 +1394,24 @@ bool GDALExtendedDataType::CopyValue(const void* pSrc,
     {
         const char* srcStrPtr;
         memcpy(&srcStrPtr, pSrc, sizeof(const char*));
-        const double dfVal = srcStrPtr == nullptr ? 0 : CPLAtof(srcStrPtr);
-        GDALCopyWords( &dfVal, GDT_Float64, 0,
-                       pDst, dstType.GetNumericDataType(), 0,
-                       1 );
+        if( dstType.GetNumericDataType() == GDT_Int64 )
+        {
+            *(static_cast<int64_t*>(pDst)) = static_cast<int64_t>(CPLAtoGIntBig(srcStrPtr));
+        }
+#if HAVE_STRTOULL
+        else if( dstType.GetNumericDataType() == GDT_UInt64 )
+        {
+            *(static_cast<uint64_t*>(pDst)) = static_cast<uint64_t>(strtoull(srcStrPtr, nullptr, 10));
+        }
+#endif
+        else
+        {
+            // FIXME GDT_UInt64
+            const double dfVal = srcStrPtr == nullptr ? 0 : CPLAtof(srcStrPtr);
+            GDALCopyWords( &dfVal, GDT_Float64, 0,
+                           pDst, dstType.GetNumericDataType(), 0,
+                           1 );
+        }
         return true;
     }
     if( srcType.GetClass() == GEDTC_COMPOUND &&
@@ -5610,6 +5635,24 @@ bool GDALMDArrayMask::IRead(const GUInt64* arrayStartIdx,
                                 bHasValidMax, dfValidMax);
             break;
 
+        case GDT_UInt64:
+            ReadInternal<std::uint64_t>(count, bufferStride, bufferDataType, pDstBuffer,
+                                pTempBuffer, oTmpBufferDT, tmpBufferStrideVector,
+                                bHasMissingValue, dfMissingValue,
+                                bHasFillValue, dfFillValue,
+                                bHasValidMin, dfValidMin,
+                                bHasValidMax, dfValidMax);
+            break;
+
+        case GDT_Int64:
+            ReadInternal<std::int64_t>(count, bufferStride, bufferDataType, pDstBuffer,
+                                pTempBuffer, oTmpBufferDT, tmpBufferStrideVector,
+                                bHasMissingValue, dfMissingValue,
+                                bHasFillValue, dfFillValue,
+                                bHasValidMin, dfValidMin,
+                                bHasValidMax, dfValidMax);
+            break;
+
         case GDT_Float32:
             ReadInternal<float>(count, bufferStride, bufferDataType, pDstBuffer,
                                 pTempBuffer, oTmpBufferDT, tmpBufferStrideVector,
diff --git a/gcore/gdalnodatamaskband.cpp b/gcore/gdalnodatamaskband.cpp
index 3a80905efa4f..004f3dcd4c39 100644
--- a/gcore/gdalnodatamaskband.cpp
+++ b/gcore/gdalnodatamaskband.cpp
@@ -102,6 +102,12 @@ static GDALDataType GetWorkDataType(GDALDataType eDataType)
         eWrkDT = GDT_Float64;
         break;
 
+      case GDT_Int64:
+      case GDT_UInt64:
+        // FIXME when we can request the nodata value as something != double
+        eWrkDT = GDT_Float64;
+        break;
+
       default:
         CPLAssert( false );
         eWrkDT = GDT_Float64;
diff --git a/gcore/gdalrasterband.cpp b/gcore/gdalrasterband.cpp
index 15574b0832a5..16dc72350260 100644
--- a/gcore/gdalrasterband.cpp
+++ b/gcore/gdalrasterband.cpp
@@ -4784,6 +4784,12 @@ inline double GetPixelValue( GDALDataType eDataType,
         case GDT_Int32:
             dfValue = static_cast<const GInt32 *>(pData)[iOffset];
             break;
+        case GDT_UInt64:
+            dfValue = static_cast<double>(static_cast<const std::uint64_t *>(pData)[iOffset]);
+            break;
+        case GDT_Int64:
+            dfValue = static_cast<double>(static_cast<const std::int64_t *>(pData)[iOffset]);
+            break;
         case GDT_Float32:
         {
             const float fValue = static_cast<const float *>(pData)[iOffset];
diff --git a/gcore/rasterio.cpp b/gcore/rasterio.cpp
index 5699c307e012..60ec8df08cc6 100644
--- a/gcore/rasterio.cpp
+++ b/gcore/rasterio.cpp
@@ -2857,6 +2857,16 @@ inline void GDALCopyWordsFromT( const T* const CPL_RESTRICT pSrcData,
                         static_cast<int*>(pDstData), nDstPixelStride,
                         nWordCount );
         break;
+    case GDT_UInt64:
+        GDALCopyWordsT( pSrcData, nSrcPixelStride,
+                        static_cast<std::uint64_t*>(pDstData), nDstPixelStride,
+                        nWordCount );
+        break;
+    case GDT_Int64:
+        GDALCopyWordsT( pSrcData, nSrcPixelStride,
+                        static_cast<std::int64_t*>(pDstData), nDstPixelStride,
+                        nWordCount );
+        break;
     case GDT_Float32:
         GDALCopyWordsT( pSrcData, nSrcPixelStride,
                         static_cast<float*>(pDstData), nDstPixelStride,
@@ -3037,6 +3047,8 @@ void GDALReplicateWord( const void * CPL_RESTRICT pSrcData,
         CASE_DUPLICATE_SIMPLE(GDT_Int16,  GInt16)
         CASE_DUPLICATE_SIMPLE(GDT_UInt32, GUInt32)
         CASE_DUPLICATE_SIMPLE(GDT_Int32,  GInt32)
+        CASE_DUPLICATE_SIMPLE(GDT_UInt64,  std::uint64_t)
+        CASE_DUPLICATE_SIMPLE(GDT_Int64,   std::int64_t)
         CASE_DUPLICATE_SIMPLE(GDT_Float32, float)
         CASE_DUPLICATE_SIMPLE(GDT_Float64, double)
 
@@ -3532,6 +3544,18 @@ GDALCopyWords64( const void * CPL_RESTRICT pSrcData,
             pDstData, eDstType, nDstPixelStride,
             nWordCount );
         break;
+    case GDT_UInt64:
+        GDALCopyWordsFromT<std::uint64_t>(
+            static_cast<const std::uint64_t *>(pSrcData), nSrcPixelStride, false,
+            pDstData, eDstType, nDstPixelStride,
+            nWordCount );
+        break;
+    case GDT_Int64:
+        GDALCopyWordsFromT<std::int64_t>(
+            static_cast<const std::int64_t *>(pSrcData), nSrcPixelStride, false,
+            pDstData, eDstType, nDstPixelStride,
+            nWordCount );
+        break;
     case GDT_Float32:
         GDALCopyWordsFromT<float>(
             static_cast<const float *>(pSrcData), nSrcPixelStride, false,
@@ -5313,6 +5337,23 @@ bool GDALBufferHasOnlyNoData( const void* pBuffer,
                                   static_cast<int32_t>(dfNoDataValue)),
                               nWidth, nHeight, nLineStride, nComponents);
     }
+    if( nBitsPerSample == 64 && nSampleFormat == GSF_UNSIGNED_INT )
+    {
+        return GDALIsValueInRange<uint64_t>(dfNoDataValue) &&
+               HasOnlyNoDataT(static_cast<const uint64_t*>(pBuffer),
+                              static_cast<uint64_t>(dfNoDataValue),
+                              nWidth, nHeight, nLineStride, nComponents);
+    }
+    if( nBitsPerSample == 64 && nSampleFormat == GSF_SIGNED_INT )
+    {
+        // Use unsigned implementation by converting the nodatavalue to
+        // unsigned
+        return GDALIsValueInRange<int64_t>(dfNoDataValue) &&
+               HasOnlyNoDataT(static_cast<const uint64_t*>(pBuffer),
+                              static_cast<uint64_t>(
+                                  static_cast<int64_t>(dfNoDataValue)),
+                              nWidth, nHeight, nLineStride, nComponents);
+    }
     if( nBitsPerSample == 32 && nSampleFormat == GSF_FLOATING_POINT )
     {
         return (std::isnan(dfNoDataValue) ||
diff --git a/swig/include/gdal.i b/swig/include/gdal.i
index 8eb72a44c826..6f5845e46f97 100644
--- a/swig/include/gdal.i
+++ b/swig/include/gdal.i
@@ -149,6 +149,8 @@ typedef enum {
     /*! Sixteen bit signed integer */           GDT_Int16 = 3,
     /*! Thirty two bit unsigned integer */      GDT_UInt32 = 4,
     /*! Thirty two bit signed integer */        GDT_Int32 = 5,
+    /*! 64 bit unsigned integer */              GDT_UInt64 = 12,
+    /*! 64 bit signed integer */                GDT_Int64 = 13,
     /*! Thirty two bit floating point */        GDT_Float32 = 6,
     /*! Sixty four bit floating point */        GDT_Float64 = 7,
     /*! Complex Int16 */                        GDT_CInt16 = 8,
diff --git a/swig/include/gdal_array.i b/swig/include/gdal_array.i
index 82496acf3aad..2635d040bc31 100644
--- a/swig/include/gdal_array.i
+++ b/swig/include/gdal_array.i
@@ -401,18 +401,22 @@ static GDALDataType NumpyTypeToGDALType(PyArrayObject *psArray)
       case NPY_FLOAT:
         return GDT_Float32;
 
-      case NPY_INT:
-      case NPY_LONG:
+      case NPY_INT32:
         return GDT_Int32;
 
-      case NPY_UINT:
-      case NPY_ULONG:
+      case NPY_UINT32:
         return GDT_UInt32;
 
-      case NPY_SHORT:
+      case NPY_INT64:
+        return GDT_Int64;
+
+      case NPY_UINT64:
+        return GDT_UInt64;
+
+      case NPY_INT16:
         return GDT_Int16;
 
-      case NPY_USHORT:
+      case NPY_UINT16:
         return GDT_UInt16;
 
       case NPY_BYTE:
@@ -1007,6 +1011,8 @@ static bool CheckNumericDataType(GDALExtendedDataTypeHS* dt)
         case GDT_UInt16: numpytype = NPY_UINT16; break;
         case GDT_Int32: numpytype = NPY_INT32; break;
         case GDT_UInt32: numpytype = NPY_UINT32; break;
+        case GDT_Int64: numpytype = NPY_INT64; break;
+        case GDT_UInt64: numpytype = NPY_UINT64; break;
         case GDT_Float32: numpytype = NPY_FLOAT32; break;
         case GDT_Float64: numpytype = NPY_FLOAT64; break;
         //case GDT_CInt16: numpytype = NPY_INT16; break;
@@ -1344,6 +1350,8 @@ codes = {gdalconst.GDT_Byte: numpy.uint8,
          gdalconst.GDT_Int16: numpy.int16,
          gdalconst.GDT_UInt32: numpy.uint32,
          gdalconst.GDT_Int32: numpy.int32,
+         gdalconst.GDT_UInt64: numpy.uint64,
+         gdalconst.GDT_Int64: numpy.int64,
          gdalconst.GDT_Float32: numpy.float32,
          gdalconst.GDT_Float64: numpy.float64,
          gdalconst.GDT_CInt16: numpy.complex64,
diff --git a/swig/include/gdalconst.i b/swig/include/gdalconst.i
index 04d892ea5b79..952f2ceb6924 100644
--- a/swig/include/gdalconst.i
+++ b/swig/include/gdalconst.i
@@ -58,6 +58,8 @@
 %constant GDT_Int16     = GDT_Int16;
 %constant GDT_UInt32    = GDT_UInt32;
 %constant GDT_Int32     = GDT_Int32;
+%constant GDT_UInt64    = GDT_UInt64;
+%constant GDT_Int64     = GDT_Int64;
 %constant GDT_Float32   = GDT_Float32;
 %constant GDT_Float64   = GDT_Float64;
 %constant GDT_CInt16    = GDT_CInt16;
diff --git a/swig/include/perl/gdal_perl.i b/swig/include/perl/gdal_perl.i
index 59bd9807fc4e..27ce290dcb36 100644
--- a/swig/include/perl/gdal_perl.i
+++ b/swig/include/perl/gdal_perl.i
@@ -481,6 +481,9 @@ sub DataTypeValueRange {
     return (-2147483648,2147483647) if $t =~/Int32/;
     return (-4294967295.0,4294967295.0) if $t =~/Float32/;
     return (-4294967295.0,4294967295.0) if $t =~/Float64/;
+    return (-4294967295.0,4294967295.0) if $t =~/Float32/;
+    return (-9223372036854775808,9223372036854775807) if $t =~/Int64/;
+    return (0,18446744073709551615) if $t =~/UInt64/;
 }
 
 sub DataTypeIsComplex {
diff --git a/swig/include/python/gdal_python.i b/swig/include/python/gdal_python.i
index 31e4b71e5b22..ef316e6fd34f 100644
--- a/swig/include/python/gdal_python.i
+++ b/swig/include/python/gdal_python.i
@@ -31,6 +31,8 @@ static int getAlignment(GDALDataType ntype)
         case GDT_Float32:
             return 4;
         case GDT_Float64:
+        case GDT_Int64:
+        case GDT_UInt64:
             return 8;
         case GDT_CInt16:
             return 2;
diff --git a/swig/python/extensions/gdal_array_wrap.cpp b/swig/python/extensions/gdal_array_wrap.cpp
index f1f44ca1f795..5287c4634359 100644
--- a/swig/python/extensions/gdal_array_wrap.cpp
+++ b/swig/python/extensions/gdal_array_wrap.cpp
@@ -3546,18 +3546,22 @@ static GDALDataType NumpyTypeToGDALType(PyArrayObject *psArray)
       case NPY_FLOAT:
         return GDT_Float32;
 
-      case NPY_INT:
-      case NPY_LONG:
+      case NPY_INT32:
         return GDT_Int32;
 
-      case NPY_UINT:
-      case NPY_ULONG:
+      case NPY_UINT32:
         return GDT_UInt32;
 
-      case NPY_SHORT:
+      case NPY_INT64:
+        return GDT_Int64;
+
+      case NPY_UINT64:
+        return GDT_UInt64;
+
+      case NPY_INT16:
         return GDT_Int16;
 
-      case NPY_USHORT:
+      case NPY_UINT16:
         return GDT_UInt16;
 
       case NPY_BYTE:
@@ -5619,6 +5623,8 @@ SWIGINTERN PyObject *_wrap_VirtualMemGetArray(PyObject *SWIGUNUSEDPARM(self), Py
       case GDT_UInt16: numpytype = NPY_UINT16; break;
       case GDT_Int32: numpytype = NPY_INT32; break;
       case GDT_UInt32: numpytype = NPY_UINT32; break;
+      case GDT_Int64: numpytype = NPY_INT64; break;
+      case GDT_UInt64: numpytype = NPY_UINT64; break;
       case GDT_Float32: numpytype = NPY_FLOAT32; break;
       case GDT_Float64: numpytype = NPY_FLOAT64; break;
       //case GDT_CInt16: numpytype = NPY_INT16; break;
diff --git a/swig/python/extensions/gdal_wrap.cpp b/swig/python/extensions/gdal_wrap.cpp
index c9164f0f1ff1..51ac44a806cd 100644
--- a/swig/python/extensions/gdal_wrap.cpp
+++ b/swig/python/extensions/gdal_wrap.cpp
@@ -2947,6 +2947,8 @@ static int getAlignment(GDALDataType ntype)
         case GDT_Float32:
             return 4;
         case GDT_Float64:
+        case GDT_Int64:
+        case GDT_UInt64:
             return 8;
         case GDT_CInt16:
             return 2;
diff --git a/swig/python/extensions/gdalconst_wrap.c b/swig/python/extensions/gdalconst_wrap.c
index 7be825a4e3d2..dae4817dd8ec 100644
--- a/swig/python/extensions/gdalconst_wrap.c
+++ b/swig/python/extensions/gdalconst_wrap.c
@@ -3523,6 +3523,8 @@ SWIG_init(void) {
   SWIG_Python_SetConstant(d, "GDT_Int16",SWIG_From_int((int)(GDT_Int16)));
   SWIG_Python_SetConstant(d, "GDT_UInt32",SWIG_From_int((int)(GDT_UInt32)));
   SWIG_Python_SetConstant(d, "GDT_Int32",SWIG_From_int((int)(GDT_Int32)));
+  SWIG_Python_SetConstant(d, "GDT_UInt64",SWIG_From_int((int)(GDT_UInt64)));
+  SWIG_Python_SetConstant(d, "GDT_Int64",SWIG_From_int((int)(GDT_Int64)));
   SWIG_Python_SetConstant(d, "GDT_Float32",SWIG_From_int((int)(GDT_Float32)));
   SWIG_Python_SetConstant(d, "GDT_Float64",SWIG_From_int((int)(GDT_Float64)));
   SWIG_Python_SetConstant(d, "GDT_CInt16",SWIG_From_int((int)(GDT_CInt16)));
diff --git a/swig/python/gdal-utils/osgeo_utils/gdal_polygonize.py b/swig/python/gdal-utils/osgeo_utils/gdal_polygonize.py
index ff9616fe0d24..883e591fc0b6 100644
--- a/swig/python/gdal-utils/osgeo_utils/gdal_polygonize.py
+++ b/swig/python/gdal-utils/osgeo_utils/gdal_polygonize.py
@@ -128,7 +128,11 @@ def gdal_polygonize(src_filename: Optional[str] = None, band_number: Union[int,
         if dst_fieldname is None:
             dst_fieldname = 'DN'
 
-        fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
+        data_type = ogr.OFTInteger
+        if srcband.DataType == gdal.GDT_Int64 or srcband.DataType == gdal.GDT_UInt64:
+            data_type = ogr.OFTInteger64
+
+        fd = ogr.FieldDefn(dst_fieldname, data_type)
         dst_layer.CreateField(fd)
         dst_field = 0
     else:
diff --git a/swig/python/osgeo/gdal_array.py b/swig/python/osgeo/gdal_array.py
index 189f71ed872c..757b93973566 100644
--- a/swig/python/osgeo/gdal_array.py
+++ b/swig/python/osgeo/gdal_array.py
@@ -141,6 +141,8 @@ def RATValuesIONumPyRead(poRAT, nField, nStart, nLength):
          gdalconst.GDT_Int16: numpy.int16,
          gdalconst.GDT_UInt32: numpy.uint32,
          gdalconst.GDT_Int32: numpy.int32,
+         gdalconst.GDT_UInt64: numpy.uint64,
+         gdalconst.GDT_Int64: numpy.int64,
          gdalconst.GDT_Float32: numpy.float32,
          gdalconst.GDT_Float64: numpy.float64,
          gdalconst.GDT_CInt16: numpy.complex64,
diff --git a/swig/python/osgeo/gdalconst.py b/swig/python/osgeo/gdalconst.py
index 5b0f10e00c42..9a43bca6e441 100644
--- a/swig/python/osgeo/gdalconst.py
+++ b/swig/python/osgeo/gdalconst.py
@@ -67,6 +67,8 @@ class _SwigNonDynamicMeta(type):
 GDT_Int16 = _gdalconst.GDT_Int16
 GDT_UInt32 = _gdalconst.GDT_UInt32
 GDT_Int32 = _gdalconst.GDT_Int32
+GDT_UInt64 = _gdalconst.GDT_UInt64
+GDT_Int64 = _gdalconst.GDT_Int64
 GDT_Float32 = _gdalconst.GDT_Float32
 GDT_Float64 = _gdalconst.GDT_Float64
 GDT_CInt16 = _gdalconst.GDT_CInt16