UVRaster: fix reprojection issue of map extent to layer projection (f…

…ixes #5501)

mapdraw.c

@@ -956,8 +956,100 @@ int msDrawVectorLayer(mapObj *map, layerObj *layer, imageObj *image)
   if(layer->transform == MS_TRUE) {
     searchrect = map->extent;
 #ifdef USE_PROJ
-    if((map->projection.numargs > 0) && (layer->projection.numargs > 0))
-      msProjectRect(&map->projection, &layer->projection, &searchrect); /* project the searchrect to source coords */
+    if((map->projection.numargs > 0) && (layer->projection.numargs > 0)) {
+      int bDone = MS_FALSE;
+
+      /* For UVRaster, it is important that the searchrect is not too large */
+      /* to avoid insufficient intermediate raster resolution, which could */
+      /* happen if we use the default code path, given potential reprojection */
+      /* issues when using a map extent that is not in the validity area of */
+      /* the layer projection. */
+      if( layer->connectiontype == MS_UVRASTER &&
+          !layer->projection.gt.need_geotransform &&
+          !(pj_is_latlong(map->projection.proj) &&
+            pj_is_latlong(layer->projection.proj)) ) {
+        rectObj layer_ori_extent;
+
+        if( msLayerGetExtent(layer, &layer_ori_extent) == MS_SUCCESS ) {
+          projectionObj map_proj;
+
+          double map_extent_minx = map->extent.minx;
+          double map_extent_miny = map->extent.miny;
+          double map_extent_maxx = map->extent.maxx;
+          double map_extent_maxy = map->extent.maxy;
+          rectObj layer_extent = layer_ori_extent;
+
+          /* Create a variant of map->projection without geotransform for */
+          /* conveniency */
+          msInitProjection(&map_proj);
+          msCopyProjection(&map_proj, &map->projection);
+          map_proj.gt.need_geotransform = MS_FALSE;
+          if( map->projection.gt.need_geotransform ) {
+            map_extent_minx = map->projection.gt.geotransform[0]
+                + map->projection.gt.geotransform[1] * map->extent.minx
+                + map->projection.gt.geotransform[2] * map->extent.miny;
+            map_extent_miny = map->projection.gt.geotransform[3]
+                + map->projection.gt.geotransform[4] * map->extent.minx
+                + map->projection.gt.geotransform[5] * map->extent.miny;
+            map_extent_maxx = map->projection.gt.geotransform[0]
+                + map->projection.gt.geotransform[1] * map->extent.maxx
+                + map->projection.gt.geotransform[2] * map->extent.maxy;
+            map_extent_maxy = map->projection.gt.geotransform[3]
+                + map->projection.gt.geotransform[4] * map->extent.maxx
+                + map->projection.gt.geotransform[5] * map->extent.maxy;
+          }
+
+          /* Reproject layer extent to map projection */
+          msProjectRect(&layer->projection, &map_proj, &layer_extent);
+
+          if( layer_extent.minx <= map_extent_minx &&
+              layer_extent.miny <= map_extent_miny &&
+              layer_extent.maxx >= map_extent_maxx &&
+              layer_extent.maxy >= map_extent_maxy ) {
+            /* do nothing special if area to map is inside layer extent */
+          }
+          else {
+            if( layer_extent.minx >= map_extent_minx &&
+                layer_extent.maxx <= map_extent_maxx &&
+                layer_extent.miny >= map_extent_miny &&
+                layer_extent.maxy <= map_extent_maxy ) {
+              /* if the area to map is larger than the layer extent, then */
+              /* use full layer extent and add some margin to reflect the */
+              /* proportion of the useful area over the requested bbox */
+              double extra_x =
+                (map_extent_maxx - map_extent_minx) /
+                  (layer_extent.maxx - layer_extent.minx) *
+                  (layer_ori_extent.maxx -  layer_ori_extent.minx);
+              double extra_y =
+                 (map_extent_maxy - map_extent_miny) /
+                  (layer_extent.maxy - layer_extent.miny) *
+                  (layer_ori_extent.maxy -  layer_ori_extent.miny);
+              searchrect.minx = layer_ori_extent.minx - extra_x / 2;
+              searchrect.maxx = layer_ori_extent.maxx + extra_x / 2;
+              searchrect.miny = layer_ori_extent.miny - extra_y / 2;
+              searchrect.maxy = layer_ori_extent.maxy + extra_y / 2;
+            }
+            else
+            {
+              /* otherwise clip the map extent with the reprojected layer */
+              /* extent */
+              searchrect.minx = MAX( map_extent_minx, layer_extent.minx );
+              searchrect.maxx = MIN( map_extent_maxx, layer_extent.maxx );
+              searchrect.miny = MAX( map_extent_miny, layer_extent.miny );
+              searchrect.maxy = MIN( map_extent_maxy, layer_extent.maxy );
+              /* and reproject into the layer projection */
+              msProjectRect(&map_proj, &layer->projection, &searchrect);
+            }
+            bDone = MS_TRUE;
+          }
+
+          msFreeProjection(&map_proj);
+        }
+      }
+
+      if( !bDone )
+        msProjectRect(&map->projection, &layer->projection, &searchrect); /* project the searchrect to source coords */
+    }
 #endif
   } else {
     searchrect.minx = searchrect.miny = 0;

mapproject.c

@@ -930,7 +930,7 @@ msProjectRectAsPolygon(projectionObj *in, projectionObj *out,
   /*      logic.                                                          */
   /* -------------------------------------------------------------------- */
   if( out && pj_is_latlong(out->proj) && in && !pj_is_latlong(in->proj)
-      && rect->maxx - rect->minx > 360.0 ) {
+      && rect->maxx - rect->minx > 360.0 && !out->gt.need_geotransform ) {
     rect->maxx = 180;
     rect->minx = -180;
   }

mapuvraster.c

@@ -698,6 +698,8 @@ int msUVRASTERLayerGetExtent(layerObj *layer, rectObj *extent)
   msTryBuildPath3(szPath, map->mappath, map->shapepath, layer->data);
   decrypted_path = msDecryptStringTokens( map, szPath );
 
+  GDALAllRegister();
+
   msAcquireLock( TLOCK_GDAL );
   if( decrypted_path ) {
     hDS = GDALOpen(decrypted_path, GA_ReadOnly );

msautotest/wxs/wms_uvraster_map_reprojection.map

@@ -0,0 +1,77 @@
+#
+# Test WMS
+#
+# REQUIRES: INPUT=GDAL OUTPUT=PNG SUPPORTS=WMS
+# 
+# RUN_PARMS: wms_uvraster_map_reprojection_extent_larger_than_layer.png [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=-85,-500,190,200&CRS=EPSG:4326&WIDTH=946&HEIGHT=583&LAYERS=test&STYLES=&FORMAT=image/png&DPI=98&MAP_RESOLUTION=98&FORMAT_OPTIONS=dpi:98&TRANSPARENT=TRUE" > [RESULT_DEMIME]
+# 
+# RUN_PARMS: wms_uvraster_map_reprojection_extent_smaller_than_layer.png [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=30,-170,60,-100&CRS=EPSG:4326&WIDTH=400&HEIGHT=200&LAYERS=test&STYLES=&FORMAT=image/png&TRANSPARENT=TRUE" > [RESULT_DEMIME]
+# 
+# RUN_PARMS: wms_uvraster_map_reprojection_extent_intersecting_layer.png [MAPSERV] QUERY_STRING="map=[MAPFILE]&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&BBOX=30,-200,60,-100&CRS=EPSG:4326&WIDTH=400&HEIGHT=200&LAYERS=test&STYLES=&FORMAT=image/png&TRANSPARENT=TRUE" > [RESULT_DEMIME]
+
+MAP
+ NAME test
+ EXTENT -180 -90 180 90
+ MAXSIZE 4096
+ SIZE 500 300
+ SHAPEPATH ./data
+
+ SYMBOL
+  NAME "arrow_wind"
+  TYPE vector
+  FILLED true
+  POINTS
+   0 1.2
+   6 1.2
+   6 3
+   10 0
+   6 -3
+   6 -1.2
+   0 -1.2
+   0 1.2
+  END
+ END
+
+ PROJECTION
+  "init=epsg:4326"
+ END
+
+ WEB
+  METADATA
+    "ows_enable_request" "*"
+  END
+ END
+
+ LAYER
+  NAME "test"
+  TYPE POINT
+  CONNECTIONTYPE uvraster
+  PROCESSING "BANDS=1,2"
+  PROCESSING "UV_SPACING=20"
+
+  PROJECTION
+   "proj=stere"
+   "lat_0=90"
+   "lat_ts=60"
+   "lon_0=249"
+   "k=90"
+   "x_0=0"
+   "y_0=0"
+   "a=6371229"
+   "b=6371229"
+   "units=m"
+   "no_defs"
+  END
+
+  DATA 'wms_uvraster_map_reprojection.tif'
+
+  CLASS
+   STYLE
+    COLOR 0 0 127
+    SYMBOL "arrow_wind"
+    ANGLE [uv_angle]
+    SIZE 12
+   END
+  END
+ END
+END