GetFeatureInfo: fix for layers using a lon_wrap PROJECTION and requested in the wrapped part of the raster

Author: rouault

maprasterquery.c

@@ -415,6 +415,9 @@ msRasterQueryByRectLow(mapObj *map, layerObj *layer, GDALDatasetH hDS,
   CPLErr      eErr;
   rasterLayerInfo *rlinfo;
   rectObj     searchrect;
+#if PROJ_VERSION_MAJOR < 6
+  int         mayNeedLonWrapAdjustment = MS_FALSE;
+#endif
 
   rlinfo = (rasterLayerInfo *) layer->layerinfo;
 
@@ -551,6 +554,16 @@ msRasterQueryByRectLow(mapObj *map, layerObj *layer, GDALDatasetH hDS,
     + sqrt( rlinfo->range_dist );
   dfAdjustedRange = dfAdjustedRange * dfAdjustedRange;
 
+#if PROJ_VERSION_MAJOR < 6
+    if( layer->project &&
+        msProjIsGeographicCRS(&(layer->projection)) &&
+        msProjIsGeographicCRS(&(map->projection)) )
+    {
+        double dfLonWrap = 0;
+        mayNeedLonWrapAdjustment = msProjectHasLonWrap(&(layer->projection), &dfLonWrap);
+    }
+#endif
+
   /* -------------------------------------------------------------------- */
   /*      Loop over all pixels determining which are "in".                */
   /* -------------------------------------------------------------------- */
@@ -576,8 +589,23 @@ msRasterQueryByRectLow(mapObj *map, layerObj *layer, GDALDatasetH hDS,
       /* coordinates if we have a hit */
       sReprojectedPixelLocation = sPixelLocation;
       if( layer->project )
+      {
+#if PROJ_VERSION_MAJOR < 6
+        /* Works around a bug in PROJ < 6 when reprojecting from a lon_wrap */
+        /* geogCRS to a geogCRS, and the input abs(longitude) is > 180. Then */
+        /* lon_wrap was ignored and the output longitude remained as the source */
+        if( mayNeedLonWrapAdjustment )
+        {
+            if( rlinfo->target_point.x < sReprojectedPixelLocation.x - 180 )
+                sReprojectedPixelLocation.x -= 360;
+            else if( rlinfo->target_point.x > sReprojectedPixelLocation.x + 180 )
+                sReprojectedPixelLocation.x += 360;
+        }
+#endif
+
         msProjectPoint( &(layer->projection), &(map->projection),
                         &sReprojectedPixelLocation);
+      }
 
       /* If we are doing QueryByShape, check against the shape now */
       if( rlinfo->searchshape != NULL ) {

msautotest/wxs/expected/wms_getfeatureinfo_lon_wrap_180_left_part_raster.xml

@@ -0,0 +1,11 @@
+GetFeatureInfo results:
+
+Layer 'test'
+  Feature 0: 
+    x = '49.4'
+    y = '-4.5'
+    value_0 = '127'
+    value_list = '127'
+    red = '127'
+    green = '127'
+    blue = '127'

msautotest/wxs/expected/wms_getfeatureinfo_lon_wrap_180_right_part_raster.xml

@@ -0,0 +1,11 @@
+GetFeatureInfo results:
+
+Layer 'test'
+  Feature 0: 
+    x = '-40.6'
+    y = '-4.5'
+    value_0 = '63'
+    value_list = '63'
+    red = '63'
+    green = '63'
+    blue = '63'

msautotest/wxs/wms_getfeatureinfo_lon_wrap_180.map

@@ -0,0 +1,59 @@
+#
+# Test WMS GetfeatureInfo with a layer with lon_wrap=180
+#
+# REQUIRES: SUPPORTS=WMS
+# 
+# RUN_PARMS: wms_getfeatureinfo_lon_wrap_180_right_part_raster.xml [MAPSERV] "QUERY_STRING=map=[MAPFILE]&service=WMS&request=GetFeatureInfo&version=1.3.0&CRS=EPSG:4326&width=200&height=200&layers=test&bbox=-90,-90,90,90&format=image/png&query_layers=test&i=50&j=100" > [RESULT_DEMIME]
+# RUN_PARMS: wms_getfeatureinfo_lon_wrap_180_left_part_raster.xml [MAPSERV] "QUERY_STRING=map=[MAPFILE]&service=WMS&request=GetFeatureInfo&version=1.3.0&CRS=EPSG:4326&width=200&height=200&layers=test&bbox=-90,-90,90,90&format=image/png&query_layers=test&i=150&j=100" > [RESULT_DEMIME]
+
+MAP
+
+NAME TEST
+STATUS ON
+SIZE 40 20
+EXTENT -180 -90 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
+
+WEB
+    METADATA
+        "ows_enable_request" "*"
+        "ows_srs" "EPSG:4326"
+        "ows_title" "WMS Getfeatureinfo"
+        "wms_onlineresource" "http://localhost/cgi-bin/mapserv?map=mymap.map"
+    END
+    IMAGEPATH '/tmp/ms_tmp/'
+    IMAGEURL '/ms_tmp/'
+END
+
+LAYER
+  NAME test
+  TYPE RASTER
+  STATUS ON
+  DATA "../gdal/data/lon_wrap_180.asc"
+  METADATA
+        "wms_title" "test"
+        "ows_srs" "EPSG:4326"
+        "ows_include_items" "all"
+  END
+  PROJECTION
+    "proj=longlat"
+    "datum=WGS84"
+    "lon_wrap=180"
+    "no_defs"
+  END
+  TEMPLATE 'blank.html'
+END
+
+END # of map file