Merge pull request OSGeo#2769 from rouault/fix_vrt_rounding

Fix & improves issues related to VRT and SrcRect/DstRect
masseyr · Jul 16, 2020 · b99148c · b99148c
2 parents 0585986 + a560bf7
commit b99148c
Show file tree

Hide file tree

Showing 3 changed files with 39 additions and 7 deletions.
diff --git a/autotest/utilities/test_gdalbuildvrt_lib.py b/autotest/utilities/test_gdalbuildvrt_lib.py
@@ -101,3 +101,14 @@ def test_gdalbuildvrt_lib_ovr():
     ds = None
     gdal.GetDriverByName('VRT').Delete(tmpfilename)
 
+
+
+def test_gdalbuildvrt_lib_te_partial_overlap():
+
+    ds = gdal.BuildVRT('', '../gcore/data/byte.tif', outputBounds=[440600, 3750060, 441860, 3751260], xRes=30, yRes=60)
+    assert ds is not None
+
+    assert ds.GetRasterBand(1).Checksum() == 8454
+    xml = ds.GetMetadata('xml:VRT')[0]
+    assert '<SrcRect xOff="0" yOff="1" xSize="19" ySize="19" />' in xml
+    assert '<DstRect xOff="4" yOff="0" xSize="38" ySize="19" />' in xml
diff --git a/gdal/apps/gdalbuildvrt_lib.cpp b/gdal/apps/gdalbuildvrt_lib.cpp
@@ -123,6 +123,7 @@ static int ArgIsNumeric( const char *pszArg )
 static int  GetSrcDstWin(DatasetProperty* psDP,
                   double we_res, double ns_res,
                   double minX, double minY, double maxX, double maxY,
+                  int nTargetXSize, int nTargetYSize,
                   double* pdfSrcXOff, double* pdfSrcYOff, double* pdfSrcXSize, double* pdfSrcYSize,
                   double* pdfDstXOff, double* pdfDstYOff, double* pdfDstXSize, double* pdfDstYSize)
 {
@@ -140,8 +141,6 @@ static int  GetSrcDstWin(DatasetProperty* psDP,
     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] < minY )
          return FALSE;
 
-    *pdfSrcXSize = psDP->nRasterXSize;
-    *pdfSrcYSize = psDP->nRasterYSize;
     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] < minX )
     {
         *pdfSrcXOff = (minX - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X]) /
@@ -166,10 +165,30 @@ static int  GetSrcDstWin(DatasetProperty* psDP,
         *pdfDstYOff =
             ((maxY - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y]) / -ns_res);
     }
-    *pdfDstXSize = (psDP->nRasterXSize *
-         psDP->adfGeoTransform[GEOTRSFRM_WE_RES] / we_res);
-    *pdfDstYSize = (psDP->nRasterYSize *
-         psDP->adfGeoTransform[GEOTRSFRM_NS_RES] / ns_res);
+
+    *pdfSrcXSize = psDP->nRasterXSize;
+    *pdfSrcYSize = psDP->nRasterYSize;
+    if( *pdfSrcXOff > 0 )
+        *pdfSrcXSize -= *pdfSrcXOff;
+    if( *pdfSrcYOff > 0 )
+        *pdfSrcYSize -= *pdfSrcYOff;
+
+    const double dfSrcToDstXSize = psDP->adfGeoTransform[GEOTRSFRM_WE_RES] / we_res;
+    *pdfDstXSize = *pdfSrcXSize * dfSrcToDstXSize;
+    const double dfSrcToDstYSize = psDP->adfGeoTransform[GEOTRSFRM_NS_RES] / ns_res;
+    *pdfDstYSize = *pdfSrcYSize * dfSrcToDstYSize;
+
+    if( *pdfDstXOff + *pdfDstXSize > nTargetXSize )
+    {
+        *pdfDstXSize = nTargetXSize - *pdfDstXOff;
+        *pdfSrcXSize = *pdfDstXSize / dfSrcToDstXSize;
+    }
+
+    if( *pdfDstYOff + *pdfDstYSize > nTargetYSize )
+    {
+        *pdfDstYSize = nTargetYSize - *pdfDstYOff;
+        *pdfSrcYSize = *pdfDstYSize / dfSrcToDstYSize;
+    }
 
     return TRUE;
 }
@@ -932,6 +951,7 @@ void VRTBuilder::CreateVRTSeparate(VRTDatasetH hVRTDS)
         {
             if ( ! GetSrcDstWin(psDatasetProperties,
                         we_res, ns_res, minX, minY, maxX, maxY,
+                        nRasterXSize, nRasterYSize,
                         &dfSrcXOff, &dfSrcYOff, &dfSrcXSize, &dfSrcYSize,
                         &dfDstXOff, &dfDstYOff, &dfDstXSize, &dfDstYSize) )
                 continue;
@@ -1062,6 +1082,7 @@ void VRTBuilder::CreateVRTNonSeparate(VRTDatasetH hVRTDS)
         double dfDstYSize;
         if ( ! GetSrcDstWin(psDatasetProperties,
                         we_res, ns_res, minX, minY, maxX, maxY,
+                        nRasterXSize, nRasterYSize,
                         &dfSrcXOff, &dfSrcYOff, &dfSrcXSize, &dfSrcYSize,
                         &dfDstXOff, &dfDstYOff, &dfDstXSize, &dfDstYSize) )
             continue;

diff --git a/gdal/frmts/vrt/vrtsources.cpp b/gdal/frmts/vrt/vrtsources.cpp
@@ -215,7 +215,7 @@ void VRTSimpleSource::SetSrcMaskBand( GDALRasterBand *poNewSrcBand )
 static double RoundIfCloseToInt(double dfValue)
 {
     double dfClosestInt = floor(dfValue + 0.5);
-    return (fabs( dfValue - dfClosestInt ) < 1e-5) ? dfClosestInt : dfValue;
+    return (fabs( dfValue - dfClosestInt ) < 1e-3) ? dfClosestInt : dfValue;
 }
 
 /************************************************************************/