WCS 2.0: fix when input raster is full world in EPSG:4326

If input raster is full world in EPSG:4326 and a GetCoverage is done with
a subsettingCRS != EPSG:4326 and with no explicit with+resolutionX (or height+
resolutionY), the request currently fails.

mapproject.c

@@ -811,13 +811,32 @@ msProjectRectAsPolygon(projectionObj *in, projectionObj *out,
   double dx, dy;
 
   /* If projecting from longlat to projected */
+  /* This hack was introduced for WFS 2.0 compliance testing, but is far */
+  /* from being perfect */
   if( out && !pj_is_latlong(out->proj) && in && pj_is_latlong(in->proj) &&
       fabs(rect->minx - -180.0) < 1e-5 && fabs(rect->miny - -90.0) < 1e-5 &&
       fabs(rect->maxx - 180.0) < 1e-5 && fabs(rect->maxy - 90.0) < 1e-5) {
-    rect->minx = -1e15;
-    rect->miny = -1e15;
-    rect->maxx = 1e15;
-    rect->maxy = 1e15;
+    pointObj pointTest;
+    pointTest.x = -180;
+    pointTest.y = 85;
+    msProjectPoint(in, out, &pointTest);
+    /* Detect if we are reprojecting from EPSG:4326 to EPSG:3857 */
+    /* and if so use more plausible bounds to avoid issues with computed */
+    /* resolution for WCS */
+    if (fabs(pointTest.x - -20037508.3427892) < 1e-5 && fabs(pointTest.y - 19971868.8804086) )
+    {
+        rect->minx = -20037508.3427892;
+        rect->miny = -20037508.3427892;
+        rect->maxx = 20037508.3427892;
+        rect->maxy = 20037508.3427892;
+    }
+    else
+    {
+        rect->minx = -1e15;
+        rect->miny = -1e15;
+        rect->maxx = 1e15;
+        rect->maxy = 1e15;
+    }
     return MS_SUCCESS;
   }
 

mapwcs20.c

@@ -4107,6 +4107,10 @@ int msWCSGetCoverage20(mapObj *map, cgiRequestObj *request,
   double x_1, x_2, y_1, y_2;
   char *coverageName, *bandlist=NULL, numbands[8];
 
+
+  int widthFromComputationInImageCRS = 0;
+  int heightFromComputationInImageCRS = 0;
+
   /* number of coverage ids should be 1 */
   if (params->ids == NULL || params->ids[0] == NULL) {
     msSetError(MS_WCSERR, "Required parameter CoverageID was not supplied.",
@@ -4268,6 +4272,42 @@ this request. Check wcs/ows_enable_request settings.", "msWCSGetCoverage20()", p
     }
 
     if(msProjectionsDiffer(&imageProj, &subsetProj)) {
+#ifdef USE_PROJ
+      /* Reprojection of source raster extent of (-180,-90,180,90) to any */
+      /* projected CRS is going to exhibit strong anomalies. So instead */
+      /* do the reverse, project the subset extent to the layer CRS, and */
+      /* see how much the subset extent takes with respect to the source */
+      /* raster extent. This is only used if output width and resolutionX (or */
+      /* (height and resolutionY) are unknown. */
+      if( ((params->width == 0 && params->resolutionX == MS_WCS20_UNBOUNDED) ||
+           (params->height == 0 && params->resolutionY == MS_WCS20_UNBOUNDED)) &&
+          (pj_is_latlong(imageProj.proj) &&
+           !pj_is_latlong(subsetProj.proj) &&
+           fabs(layer->extent.minx - -180.0) < 1e-5 &&
+           fabs(layer->extent.miny - -90.0) < 1e-5 &&
+           fabs(layer->extent.maxx - 180.0) < 1e-5 &&
+           fabs(layer->extent.maxy - 90.0) < 1e-5) )
+      {
+          rectObj subsetInImageProj = subsets;
+          if( msProjectRect(&subsetProj, &imageProj, &(subsetInImageProj)) == MS_SUCCESS )
+          {
+            subsetInImageProj.minx = MS_MAX(subsetInImageProj.minx, layer->extent.minx);
+            subsetInImageProj.miny = MS_MAX(subsetInImageProj.miny, layer->extent.miny);
+            subsetInImageProj.maxx = MS_MIN(subsetInImageProj.maxx, layer->extent.maxx);
+            subsetInImageProj.maxy = MS_MIN(subsetInImageProj.maxy, layer->extent.maxy);
+            {
+                double total = ABS(layer->extent.maxx - layer->extent.minx);
+                double part = ABS(subsetInImageProj.maxx - subsetInImageProj.minx);
+                widthFromComputationInImageCRS = MS_NINT((part * map->width) / total);
+            }
+            {
+                double total = ABS(layer->extent.maxy - layer->extent.miny);
+                double part = ABS(subsetInImageProj.maxy - subsetInImageProj.miny);
+                heightFromComputationInImageCRS = MS_NINT((part * map->height) / total);
+            }
+          }
+      }
+#endif
       msProjectRect(&imageProj, &subsetProj, &(layer->extent));
       map->extent = layer->extent;
       msFreeProjection(&(map->projection));
@@ -4318,7 +4358,9 @@ this request. Check wcs/ows_enable_request settings.", "msWCSGetCoverage20()", p
   } else if(params->resolutionX != MS_WCS20_UNBOUNDED) {
     params->width = MS_NINT((bbox.maxx - bbox.minx) / params->resolutionX);
   } else {
-    if(ABS(bbox.maxx - bbox.minx) != ABS(map->extent.maxx - map->extent.minx)) {
+    if( widthFromComputationInImageCRS != 0 ) {
+      params->width = widthFromComputationInImageCRS;
+    } else if(ABS(bbox.maxx - bbox.minx) != ABS(map->extent.maxx - map->extent.minx)) {
       double total = ABS(map->extent.maxx - map->extent.minx),
              part = ABS(bbox.maxx - bbox.minx);
       params->width = MS_NINT((part * map->width) / total);
@@ -4340,7 +4382,9 @@ this request. Check wcs/ows_enable_request settings.", "msWCSGetCoverage20()", p
   } else if(params->resolutionY != MS_WCS20_UNBOUNDED) {
     params->height = MS_NINT((bbox.maxy - bbox.miny) / params->resolutionY);
   } else {
-    if(ABS(bbox.maxy - bbox.miny) != ABS(map->extent.maxy - map->extent.miny)) {
+    if( heightFromComputationInImageCRS != 0 ) {
+      params->height = heightFromComputationInImageCRS;
+    } else if(ABS(bbox.maxy - bbox.miny) != ABS(map->extent.maxy - map->extent.miny)) {
       double total = ABS(map->extent.maxy - map->extent.miny),
              part = ABS(bbox.maxy - bbox.miny);
       params->height = MS_NINT((part * map->height) / total);