diff --git a/alg/gdalwarper.cpp b/alg/gdalwarper.cpp
index 079211ae4750..ca78650897b3 100644
--- a/alg/gdalwarper.cpp
+++ b/alg/gdalwarper.cpp
@@ -1536,6 +1536,32 @@ void CPL_STDCALL GDALWarpResolveWorkingDataType(GDALWarpOptions *psOptions)
 
     psOptions->eWorkingDataType = GDT_Byte;
 
+    // If none of the provided input nodata values can be represented in the
+    // data type of the corresponding source band, ignore them.
+    if (psOptions->hSrcDS && psOptions->padfSrcNoDataReal)
+    {
+        int nCountInvalidSrcNoDataReal = 0;
+        for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
+        {
+            GDALRasterBandH hSrcBand = GDALGetRasterBand(
+                psOptions->hSrcDS, psOptions->panSrcBands[iBand]);
+
+            if (hSrcBand &&
+                !GDALIsValueExactAs(psOptions->padfSrcNoDataReal[iBand],
+                                    GDALGetRasterDataType(hSrcBand)))
+            {
+                nCountInvalidSrcNoDataReal++;
+            }
+        }
+        if (nCountInvalidSrcNoDataReal == psOptions->nBandCount)
+        {
+            CPLFree(psOptions->padfSrcNoDataReal);
+            psOptions->padfSrcNoDataReal = nullptr;
+            CPLFree(psOptions->padfSrcNoDataImag);
+            psOptions->padfSrcNoDataImag = nullptr;
+        }
+    }
+
     for (int iBand = 0; iBand < psOptions->nBandCount; iBand++)
     {
         if (psOptions->hDstDS != nullptr)
diff --git a/autotest/cpp/test_alg.cpp b/autotest/cpp/test_alg.cpp
index d8c5d64a4b44..6d58a164164e 100644
--- a/autotest/cpp/test_alg.cpp
+++ b/autotest/cpp/test_alg.cpp
@@ -151,6 +151,35 @@ TEST_F(test_alg, GDALWarpResolveWorkingDataType_padfSrcNoDataReal)
     GDALDestroyWarpOptions(psOptions);
 }
 
+// GDALWarpResolveWorkingDataType: effect of padfSrcNoDataReal
+TEST_F(test_alg, GDALWarpResolveWorkingDataType_padfSrcNoDataReal_with_band)
+{
+    GDALDatasetUniquePtr poDS(GDALDriver::FromHandle(GDALGetDriverByName("MEM"))
+                                  ->Create("", 1, 1, 1, GDT_Byte, nullptr));
+    GDALWarpOptions *psOptions = GDALCreateWarpOptions();
+    psOptions->hSrcDS = GDALDataset::ToHandle(poDS.get());
+    psOptions->nBandCount = 1;
+    psOptions->panSrcBands =
+        static_cast<int *>(CPLMalloc(psOptions->nBandCount * sizeof(int)));
+    psOptions->panSrcBands[0] = 1;
+    psOptions->padfSrcNoDataReal =
+        static_cast<double *>(CPLMalloc(sizeof(double)));
+    psOptions->padfSrcNoDataReal[0] = 0.0;
+    GDALWarpResolveWorkingDataType(psOptions);
+    EXPECT_EQ(psOptions->eWorkingDataType, GDT_Byte);
+
+    psOptions->padfSrcNoDataReal[0] = -1.0;
+    GDALWarpResolveWorkingDataType(psOptions);
+    EXPECT_EQ(psOptions->eWorkingDataType, GDT_Byte);
+
+    psOptions->eWorkingDataType = GDT_Unknown;
+    psOptions->padfSrcNoDataReal[0] = 2.0;
+    GDALWarpResolveWorkingDataType(psOptions);
+    EXPECT_EQ(psOptions->eWorkingDataType, GDT_Byte);
+
+    GDALDestroyWarpOptions(psOptions);
+}
+
 // GDALWarpResolveWorkingDataType: effect of padfSrcNoDataImag
 TEST_F(test_alg, GDALWarpResolveWorkingDataType_padfSrcNoDataImag)
 {
diff --git a/autotest/cpp/test_gdal.cpp b/autotest/cpp/test_gdal.cpp
index be310432af49..c5046797e65f 100644
--- a/autotest/cpp/test_gdal.cpp
+++ b/autotest/cpp/test_gdal.cpp
@@ -260,12 +260,74 @@ TEST_F(test_gdal, GDALDataTypeUnion_special_cases)
                                false /* complex */),
               GDT_Int64);
 
-    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -128, 0), GDT_Int16);
-    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32768, 0), GDT_Int16);
-    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32769, 0), GDT_Int32);
-    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999, 0), GDT_Float32);
-    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999.9876, 0),
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -128, false), GDT_Int16);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32768, false), GDT_Int16);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Byte, -32769, false), GDT_Int32);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int8, 127, false), GDT_Int8);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int8, 128, false), GDT_Int16);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int16, 32767, false), GDT_Int16);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int16, 32768, false), GDT_Int32);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt16, 65535, false), GDT_UInt16);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt16, 65536, false), GDT_UInt32);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int32, INT32_MAX, false),
+              GDT_Int32);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Int32, INT32_MAX + 1.0, false),
+              GDT_Int64);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt32, UINT32_MAX, false),
+              GDT_UInt32);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_UInt32, UINT32_MAX + 1.0, false),
+              GDT_UInt64);
+
+    // (1 << 63) - 1024
+    EXPECT_EQ(
+        GDALDataTypeUnionWithValue(GDT_Int64, 9223372036854774784.0, false),
+        GDT_Int64);
+    // (1 << 63) - 512
+    EXPECT_EQ(
+        GDALDataTypeUnionWithValue(GDT_Int64, 9223372036854775296.0, false),
+        GDT_Float64);
+
+    // (1 << 64) - 2048
+    EXPECT_EQ(
+        GDALDataTypeUnionWithValue(GDT_UInt64, 18446744073709549568.0, false),
+        GDT_UInt64);
+    // (1 << 64) + 4096
+    EXPECT_EQ(
+        GDALDataTypeUnionWithValue(GDT_UInt64, 18446744073709555712.0, false),
+        GDT_Float64);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999, false),
+              GDT_Float32);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float32, -99999.9876, false),
+              GDT_Float64);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(
+                  GDT_Float32, std::numeric_limits<double>::quiet_NaN(), false),
+              GDT_Float32);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(
+                  GDT_Float32, -std::numeric_limits<double>::infinity(), false),
+              GDT_Float32);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(
+                  GDT_Float32, -std::numeric_limits<double>::infinity(), false),
+              GDT_Float32);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Float64, -99999.9876, false),
+              GDT_Float64);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(
+                  GDT_Float64, std::numeric_limits<double>::quiet_NaN(), false),
               GDT_Float64);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(
+                  GDT_Float64, -std::numeric_limits<double>::infinity(), false),
+              GDT_Float64);
+    EXPECT_EQ(GDALDataTypeUnionWithValue(
+                  GDT_Float64, -std::numeric_limits<double>::infinity(), false),
+              GDT_Float64);
+
+    EXPECT_EQ(GDALDataTypeUnionWithValue(GDT_Unknown, 0, false), GDT_Byte);
 }
 
 // Test GDALAdjustValueToDataType()
@@ -807,6 +869,106 @@ TEST_F(test_gdal, GDALIsValueExactAs)
         GDALIsValueExactAs<double>(std::numeric_limits<double>::quiet_NaN()));
 }
 
+// Test GDALIsValueExactAs()
+TEST_F(test_gdal, GDALIsValueExactAs_C_func)
+{
+    EXPECT_TRUE(GDALIsValueExactAs(0, GDT_Byte));
+    EXPECT_TRUE(GDALIsValueExactAs(255, GDT_Byte));
+    EXPECT_FALSE(GDALIsValueExactAs(-1, GDT_Byte));
+    EXPECT_FALSE(GDALIsValueExactAs(256, GDT_Byte));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Byte));
+
+    EXPECT_TRUE(GDALIsValueExactAs(-128, GDT_Int8));
+    EXPECT_TRUE(GDALIsValueExactAs(127, GDT_Int8));
+    EXPECT_FALSE(GDALIsValueExactAs(-129, GDT_Int8));
+    EXPECT_FALSE(GDALIsValueExactAs(128, GDT_Int8));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int8));
+
+    EXPECT_TRUE(GDALIsValueExactAs(0, GDT_UInt16));
+    EXPECT_TRUE(GDALIsValueExactAs(65535, GDT_UInt16));
+    EXPECT_FALSE(GDALIsValueExactAs(-1, GDT_UInt16));
+    EXPECT_FALSE(GDALIsValueExactAs(65536, GDT_UInt16));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt16));
+
+    EXPECT_TRUE(GDALIsValueExactAs(-32768, GDT_Int16));
+    EXPECT_TRUE(GDALIsValueExactAs(32767, GDT_Int16));
+    EXPECT_FALSE(GDALIsValueExactAs(-32769, GDT_Int16));
+    EXPECT_FALSE(GDALIsValueExactAs(32768, GDT_Int16));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int16));
+
+    EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<uint32_t>::lowest(),
+                                   GDT_UInt32));
+    EXPECT_TRUE(
+        GDALIsValueExactAs(std::numeric_limits<uint32_t>::max(), GDT_UInt32));
+    EXPECT_FALSE(GDALIsValueExactAs(
+        std::numeric_limits<uint32_t>::lowest() - 1.0, GDT_UInt32));
+    EXPECT_FALSE(GDALIsValueExactAs(std::numeric_limits<uint32_t>::max() + 1.0,
+                                    GDT_UInt32));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt32));
+
+    EXPECT_TRUE(
+        GDALIsValueExactAs(std::numeric_limits<int32_t>::lowest(), GDT_Int32));
+    EXPECT_TRUE(
+        GDALIsValueExactAs(std::numeric_limits<int32_t>::max(), GDT_Int32));
+    EXPECT_FALSE(GDALIsValueExactAs(
+        std::numeric_limits<int32_t>::lowest() - 1.0, GDT_Int32));
+    EXPECT_FALSE(GDALIsValueExactAs(std::numeric_limits<int32_t>::max() + 1.0,
+                                    GDT_Int32));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int32));
+
+    EXPECT_TRUE(GDALIsValueExactAs(
+        static_cast<double>(std::numeric_limits<uint64_t>::lowest()),
+        GDT_UInt64));
+    // (1 << 64) - 2048
+    EXPECT_TRUE(GDALIsValueExactAs(18446744073709549568.0, GDT_UInt64));
+    EXPECT_FALSE(GDALIsValueExactAs(
+        static_cast<double>(std::numeric_limits<uint64_t>::lowest()) - 1.0,
+        GDT_UInt64));
+    // (1 << 64)
+    EXPECT_FALSE(GDALIsValueExactAs(18446744073709551616.0, GDT_UInt64));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_UInt64));
+
+    EXPECT_TRUE(GDALIsValueExactAs(
+        static_cast<double>(std::numeric_limits<int64_t>::lowest()),
+        GDT_Int64));
+    // (1 << 63) - 1024
+    EXPECT_TRUE(GDALIsValueExactAs(9223372036854774784.0, GDT_Int64));
+    EXPECT_FALSE(GDALIsValueExactAs(
+        static_cast<double>(std::numeric_limits<int64_t>::lowest()) - 2048.0,
+        GDT_Int64));
+    // (1 << 63) - 512
+    EXPECT_FALSE(GDALIsValueExactAs(9223372036854775296.0, GDT_Int64));
+    EXPECT_FALSE(GDALIsValueExactAs(0.5, GDT_Int64));
+
+    EXPECT_TRUE(
+        GDALIsValueExactAs(-std::numeric_limits<float>::max(), GDT_Float32));
+    EXPECT_TRUE(
+        GDALIsValueExactAs(std::numeric_limits<float>::max(), GDT_Float32));
+    EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits<float>::infinity(),
+                                   GDT_Float32));
+    EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<float>::infinity(),
+                                   GDT_Float32));
+    EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<double>::quiet_NaN(),
+                                   GDT_Float32));
+    EXPECT_TRUE(
+        !GDALIsValueExactAs(-std::numeric_limits<double>::max(), GDT_Float32));
+    EXPECT_TRUE(
+        !GDALIsValueExactAs(std::numeric_limits<double>::max(), GDT_Float32));
+
+    EXPECT_TRUE(GDALIsValueExactAs(-std::numeric_limits<double>::infinity(),
+                                   GDT_Float64));
+    EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<double>::infinity(),
+                                   GDT_Float64));
+    EXPECT_TRUE(
+        GDALIsValueExactAs(-std::numeric_limits<double>::max(), GDT_Float64));
+    EXPECT_TRUE(
+        GDALIsValueExactAs(std::numeric_limits<double>::max(), GDT_Float64));
+    EXPECT_TRUE(GDALIsValueExactAs(std::numeric_limits<double>::quiet_NaN(),
+                                   GDT_Float64));
+
+    EXPECT_TRUE(GDALIsValueExactAs(0, GDT_CInt16));
+}
+
 #ifdef _MSC_VER
 #pragma warning(pop)
 #endif
diff --git a/autotest/gdrivers/vrtwarp.py b/autotest/gdrivers/vrtwarp.py
index 3f37e1d0a6c1..9b843aa81342 100755
--- a/autotest/gdrivers/vrtwarp.py
+++ b/autotest/gdrivers/vrtwarp.py
@@ -732,3 +732,31 @@ def test_vrtwarp_irasterio_optim_window_splitting():
     with gdaltest.config_option("GDAL_VRT_WARP_USE_DATASET_RASTERIO", "NO"):
         expected_data = warped_vrt_ds.ReadRaster()
     assert warped_vrt_ds.ReadRaster() == expected_data
+
+
+###############################################################################
+# Test gdal.AutoCreateWarpedVRT() on a Int16 band with nodata = 32767
+
+
+def test_vrtwarp_autocreatewarpedvrt_int16_nodata_32767():
+
+    ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1, gdal.GDT_Int16)
+    ds.SetGeoTransform([0, 1, 0, 0, 0, -1])
+    ds.GetRasterBand(1).SetNoDataValue(32767)
+    vrt_ds = gdal.AutoCreateWarpedVRT(ds)
+    assert vrt_ds.GetRasterBand(1).DataType == gdal.GDT_Int16
+    assert vrt_ds.GetRasterBand(1).GetNoDataValue() == 32767
+
+
+###############################################################################
+# Test gdal.AutoCreateWarpedVRT() on a source nodata value that does not fit
+# the source band type
+
+
+def test_vrtwarp_autocreatewarpedvrt_invalid_nodata():
+
+    ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1, gdal.GDT_Byte)
+    ds.SetGeoTransform([0, 1, 0, 0, 0, -1])
+    ds.GetRasterBand(1).SetNoDataValue(-9999)
+    vrt_ds = gdal.AutoCreateWarpedVRT(ds)
+    assert vrt_ds.GetRasterBand(1).DataType == gdal.GDT_Byte
diff --git a/doc/source/spelling_wordlist.txt b/doc/source/spelling_wordlist.txt
index b71562cf4c43..056606fb35bf 100644
--- a/doc/source/spelling_wordlist.txt
+++ b/doc/source/spelling_wordlist.txt
@@ -661,6 +661,7 @@ dfTargetHeight
 dfTimeout
 dfTolerance
 dfURY
+dfValue
 dfVisibleVal
 dfXOff
 dfXSize
diff --git a/frmts/vrt/vrtwarped.cpp b/frmts/vrt/vrtwarped.cpp
index b96cd5071183..ec39b7f1263c 100644
--- a/frmts/vrt/vrtwarped.cpp
+++ b/frmts/vrt/vrtwarped.cpp
@@ -180,6 +180,9 @@ GDALDatasetH CPL_STDCALL GDALAutoCreateWarpedVRTEx(
     if (psWO->padfSrcNoDataReal == nullptr &&
         psWO->padfDstNoDataReal == nullptr && psWO->nSrcAlphaBand == 0)
     {
+        // If none of the provided input nodata values can be represented in the
+        // data type of the corresponding source band, ignore them.
+        int nCountInvalidSrcNoDataReal = 0;
         for (int i = 0; i < psWO->nBandCount; i++)
         {
             GDALRasterBandH rasterBand =
@@ -189,22 +192,42 @@ GDALDatasetH CPL_STDCALL GDALAutoCreateWarpedVRTEx(
             double noDataValue =
                 GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);
 
-            if (hasNoDataValue)
+            if (hasNoDataValue &&
+                !GDALIsValueExactAs(noDataValue,
+                                    GDALGetRasterDataType(rasterBand)))
             {
-                // Check if the nodata value is out of range
-                int bClamped = FALSE;
-                int bRounded = FALSE;
-                CPL_IGNORE_RET_VAL(GDALAdjustValueToDataType(
-                    GDALGetRasterDataType(rasterBand), noDataValue, &bClamped,
-                    &bRounded));
-                if (!bClamped)
+                nCountInvalidSrcNoDataReal++;
+            }
+        }
+
+        if (nCountInvalidSrcNoDataReal != psWO->nBandCount)
+        {
+            for (int i = 0; i < psWO->nBandCount; i++)
+            {
+                GDALRasterBandH rasterBand =
+                    GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[i]);
+
+                int hasNoDataValue;
+                double noDataValue =
+                    GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);
+
+                if (hasNoDataValue)
                 {
-                    GDALWarpInitNoDataReal(psWO, -1e10);
-                    if (psWO->padfSrcNoDataReal != nullptr &&
-                        psWO->padfDstNoDataReal != nullptr)
+                    // Check if the nodata value is out of range
+                    int bClamped = FALSE;
+                    int bRounded = FALSE;
+                    CPL_IGNORE_RET_VAL(GDALAdjustValueToDataType(
+                        GDALGetRasterDataType(rasterBand), noDataValue,
+                        &bClamped, &bRounded));
+                    if (!bClamped)
                     {
-                        psWO->padfSrcNoDataReal[i] = noDataValue;
-                        psWO->padfDstNoDataReal[i] = noDataValue;
+                        GDALWarpInitNoDataReal(psWO, -1e10);
+                        if (psWO->padfSrcNoDataReal != nullptr &&
+                            psWO->padfDstNoDataReal != nullptr)
+                        {
+                            psWO->padfSrcNoDataReal[i] = noDataValue;
+                            psWO->padfDstNoDataReal[i] = noDataValue;
+                        }
                     }
                 }
             }
diff --git a/gcore/gdal.h b/gcore/gdal.h
index 21d1a4da18d7..565133372c0c 100644
--- a/gcore/gdal.h
+++ b/gcore/gdal.h
@@ -100,6 +100,7 @@ GDALDataType CPL_DLL CPL_STDCALL GDALFindDataTypeForValue(double dValue,
                                                           int bComplex);
 double CPL_DLL GDALAdjustValueToDataType(GDALDataType eDT, double dfValue,
                                          int *pbClamped, int *pbRounded);
+bool CPL_DLL GDALIsValueExactAs(double dfValue, GDALDataType eDT);
 GDALDataType CPL_DLL CPL_STDCALL GDALGetNonComplexDataType(GDALDataType);
 int CPL_DLL CPL_STDCALL GDALDataTypeIsConversionLossy(GDALDataType eTypeFrom,
                                                       GDALDataType eTypeTo);
diff --git a/gcore/gdal_misc.cpp b/gcore/gdal_misc.cpp
index af3d01d6e96b..8daa761b33de 100644
--- a/gcore/gdal_misc.cpp
+++ b/gcore/gdal_misc.cpp
@@ -140,6 +140,11 @@ GDALDataType CPL_STDCALL GDALDataTypeUnion(GDALDataType eType1,
                                            GDALDataType eType2)
 
 {
+    if (eType1 == GDT_Unknown)
+        return eType2;
+    if (eType2 == GDT_Unknown)
+        return eType1;
+
     const int panBits[] = {GetDataTypeElementSizeBits(eType1),
                            GetDataTypeElementSizeBits(eType2)};
 
@@ -169,19 +174,89 @@ GDALDataType CPL_STDCALL GDALDataTypeUnion(GDALDataType eType1,
  * \brief Union a data type with the one found for a value
  *
  * @param eDT the first data type
- * @param dValue the value for which to find a data type and union with eDT
+ * @param dfValue the value for which to find a data type and union with eDT
  * @param bComplex if the value is complex
  *
- * @return a data type able to express eDT and dValue.
+ * @return a data type able to express eDT and dfValue.
  * @since GDAL 2.3
  */
 GDALDataType CPL_STDCALL GDALDataTypeUnionWithValue(GDALDataType eDT,
-                                                    double dValue, int bComplex)
+                                                    double dfValue,
+                                                    int bComplex)
 {
-    if (eDT == GDT_Float32 && !bComplex && static_cast<float>(dValue) == dValue)
-        return eDT;
+    if (!bComplex && !GDALDataTypeIsComplex(eDT))
+    {
+        switch (eDT)
+        {
+            case GDT_Byte:
+            {
+                if (GDALIsValueExactAs<uint8_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_Int8:
+            {
+                if (GDALIsValueExactAs<int8_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_UInt16:
+            {
+                if (GDALIsValueExactAs<uint16_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_Int16:
+            {
+                if (GDALIsValueExactAs<int16_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_UInt32:
+            {
+                if (GDALIsValueExactAs<uint32_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_Int32:
+            {
+                if (GDALIsValueExactAs<int32_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_UInt64:
+            {
+                if (GDALIsValueExactAs<uint64_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_Int64:
+            {
+                if (GDALIsValueExactAs<int64_t>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_Float32:
+            {
+                if (GDALIsValueExactAs<float>(dfValue))
+                    return eDT;
+                break;
+            }
+            case GDT_Float64:
+            {
+                return eDT;
+            }
+            case GDT_Unknown:
+            case GDT_CInt16:
+            case GDT_CInt32:
+            case GDT_CFloat32:
+            case GDT_CFloat64:
+            case GDT_TypeCount:
+                break;
+        }
+    }
 
-    const GDALDataType eDT2 = GDALFindDataTypeForValue(dValue, bComplex);
+    const GDALDataType eDT2 = GDALFindDataTypeForValue(dfValue, bComplex);
     return GDALDataTypeUnion(eDT, eDT2);
 }
 
@@ -313,7 +388,12 @@ GDALDataType CPL_STDCALL GDALFindDataType(int nBits, int bSigned, int bFloating,
  */
 GDALDataType CPL_STDCALL GDALFindDataTypeForValue(double dValue, int bComplex)
 {
-    const bool bFloating = round(dValue) != dValue;
+    const bool bFloating =
+        round(dValue) != dValue ||
+        dValue >
+            static_cast<double>(std::numeric_limits<std::uint64_t>::max()) ||
+        dValue <
+            static_cast<double>(std::numeric_limits<std::int64_t>::lowest());
     const bool bSigned = bFloating || dValue < 0;
     const int nBits = GetMinBitsForValue(dValue);
 
@@ -868,6 +948,58 @@ double GDALAdjustValueToDataType(GDALDataType eDT, double dfValue,
     return dfValue;
 }
 
+/************************************************************************/
+/*                         GDALIsValueExactAs()                         */
+/************************************************************************/
+
+/**
+ * \brief Check whether the provided value can be exactly represented in a
+ * data type.
+ *
+ * Only implemented for non-complex data types
+ *
+ * @param dfValue value to check.
+ * @param eDT target data type.
+ *
+ * @return true if the provided value can be exactly represented in the
+ * data type.
+ * @since GDAL 3.10
+ */
+bool GDALIsValueExactAs(double dfValue, GDALDataType eDT)
+{
+    switch (eDT)
+    {
+        case GDT_Byte:
+            return GDALIsValueExactAs<uint8_t>(dfValue);
+        case GDT_Int8:
+            return GDALIsValueExactAs<int8_t>(dfValue);
+        case GDT_UInt16:
+            return GDALIsValueExactAs<uint16_t>(dfValue);
+        case GDT_Int16:
+            return GDALIsValueExactAs<int16_t>(dfValue);
+        case GDT_UInt32:
+            return GDALIsValueExactAs<uint32_t>(dfValue);
+        case GDT_Int32:
+            return GDALIsValueExactAs<int32_t>(dfValue);
+        case GDT_UInt64:
+            return GDALIsValueExactAs<uint64_t>(dfValue);
+        case GDT_Int64:
+            return GDALIsValueExactAs<int64_t>(dfValue);
+        case GDT_Float32:
+            return GDALIsValueExactAs<float>(dfValue);
+        case GDT_Float64:
+            return true;
+        case GDT_Unknown:
+        case GDT_CInt16:
+        case GDT_CInt32:
+        case GDT_CFloat32:
+        case GDT_CFloat64:
+        case GDT_TypeCount:
+            break;
+    }
+    return true;
+}
+
 /************************************************************************/
 /*                        GDALGetNonComplexDataType()                */
 /************************************************************************/