msTransformMapToSource(): extra fixes to handle lon_wrap=180 and sour…

…ce dataset whose longitude range is larger than 360 deg
MapServer · Jan 3, 2020 · 421a1fcecfa3bf5fa7ee944ea41b30eeaab734f6 · 421a1fc
1 parent 72c7845
commit 421a1fcecfa3bf5fa7ee944ea41b30eeaab734f6
Show file tree

Hide file tree

Showing 5 changed files with 70 additions and 7 deletions.
diff --git a/mapresample.c b/mapresample.c
@@ -1145,6 +1145,7 @@ static int msTransformMapToSource( int nDstXSize, int nDstYSize,
           double dfYMinOut = 0.0;
           double dfXMaxOut = 0.0;
           double dfYMaxOut = 0.0;
+          const double dfHalfRes = adfDstGeoTransform[1] / 2;
 
           /* Find out average y coordinate in src projection */
           for( i = 0; i < nSamples; i++ ) {
@@ -1183,7 +1184,7 @@ static int msTransformMapToSource( int nDstXSize, int nDstYSize,
                         2, 1, x2, y2, z2 );
           msReleaseLock( TLOCK_PROJ );
 
-          if( x2[0] >= dfXMinOut && x2[0] <= dfXMaxOut &&
+          if( x2[0] >= dfXMinOut - dfHalfRes && x2[0] <= dfXMaxOut + dfHalfRes &&
               y2[0] >= dfYMinOut && y2[0] <= dfYMaxOut )
           {
                 double x_out =      adfInvSrcGeoTransform[0]
@@ -1193,8 +1194,8 @@ static int msTransformMapToSource( int nDstXSize, int nDstYSize,
                             +   (dfLonWrap-180)*adfInvSrcGeoTransform[4]
                             +   dfY*adfInvSrcGeoTransform[5];
 
-                /* Does the raster cover a whole 360 deg range ? */
-                if( nSrcXSize == (int)(adfInvSrcGeoTransform[1] * 360 + 0.5) )
+                /* Does the raster cover, at least, a whole 360 deg range ? */
+                if( nSrcXSize >= (int)(adfInvSrcGeoTransform[1] * 360) )
                 {
                     psSrcExtent->minx = 0;
                     psSrcExtent->maxx = nSrcXSize;
@@ -1208,8 +1209,8 @@ static int msTransformMapToSource( int nDstXSize, int nDstYSize,
                 psSrcExtent->maxy = MS_MAX(psSrcExtent->maxy, y_out);
           }
 
-          if( x2[1] >= dfXMinOut && x2[1] <= dfXMaxOut &&
-              x2[1] >= dfYMinOut && y2[1] <= dfYMaxOut )
+          if( x2[1] >= dfXMinOut - dfHalfRes && x2[1] <= dfXMaxOut + dfHalfRes &&
+              y2[1] >= dfYMinOut && y2[1] <= dfYMaxOut )
           {
                 double x_out =      adfInvSrcGeoTransform[0]
                             +   (dfLonWrap+180)*adfInvSrcGeoTransform[1]
@@ -1218,8 +1219,8 @@ static int msTransformMapToSource( int nDstXSize, int nDstYSize,
                             +   (dfLonWrap+180)*adfInvSrcGeoTransform[4]
                             +   dfY*adfInvSrcGeoTransform[5];
 
-                /* Does the raster cover a whole 360 deg range ? */
-                if( nSrcXSize == (int)(adfInvSrcGeoTransform[1] * 360 + 0.5) )
+                /* Does the raster cover, at least, a whole 360 deg range ? */
+                if( nSrcXSize >= (int)(adfInvSrcGeoTransform[1] * 360) )
                 {
                     psSrcExtent->minx = 0;
                     psSrcExtent->maxx = nSrcXSize;

diff --git a/msautotest/gdal/data/lon_wrap_180_covering_more_than_360deg.asc b/msautotest/gdal/data/lon_wrap_180_covering_more_than_360deg.asc
@@ -0,0 +1,25 @@
+ncols        41
+nrows        20
+xllcorner    -4.5
+yllcorner    -90.000000000000
+cellsize     9.000000000000
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
+ 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 127 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 63 127 127
diff --git a/msautotest/gdal/data/lon_wrap_180_covering_more_than_360deg.prj b/msautotest/gdal/data/lon_wrap_180_covering_more_than_360deg.prj
@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
diff --git a/msautotest/gdal/expected/lon_wrap_180_covering_more_than_360deg.png b/msautotest/gdal/expected/lon_wrap_180_covering_more_than_360deg.png
diff --git a/msautotest/gdal/lon_wrap_180_covering_more_than_360deg.map b/msautotest/gdal/lon_wrap_180_covering_more_than_360deg.map
@@ -0,0 +1,36 @@
+# Test behaviour with a projection with lon_wrap=180
+
+MAP
+
+NAME TEST
+STATUS ON
+SIZE 40 20
+EXTENT -180 -89.9999 180 90
+IMAGECOLOR 0 0 0
+
+ PROJECTION
+  "init=epsg:4326"
+ END
+
+OUTPUTFORMAT
+  NAME "png"
+  DRIVER AGG/PNG
+  MIMETYPE "image/png"
+  IMAGEMODE RGBA
+  EXTENSION "png"
+END
+
+LAYER
+  NAME test
+  TYPE RASTER
+  STATUS default
+  DATA "data/lon_wrap_180_covering_more_than_360deg.asc"
+  PROJECTION
+    "proj=longlat"
+    "datum=WGS84"
+    "lon_wrap=180"
+    "no_defs"
+  END
+END
+
+END # of map file