Merge pull request #5122 from rouault/fix_4461

WMTS: disable clipping with TileMatrixSetLimits by default ...

autotest/cpp/test_ogr.cpp

@@ -1913,4 +1913,32 @@ namespace tut
         }
     }
 
+    // Test OGREnvelope
+    template<>
+    template<>
+    void object::test<22>()
+    {
+        OGREnvelope s1;
+        ensure(!s1.IsInit());
+        {
+            OGREnvelope s2(s1);
+            ensure(s1 == s2);
+            ensure(!(s1 != s2));
+        }
+
+        s1.MinX = 0;
+        s1.MinY = 1;
+        s1.MaxX = 2;
+        s1.MaxX = 3;
+        ensure(s1.IsInit());
+        {
+            OGREnvelope s2(s1);
+            ensure(s1 == s2);
+            ensure(!(s1 != s2));
+            s2.MinX += 1;
+            ensure(s1 != s2);
+            ensure(!(s1 == s2));
+        }
+    }
+
 } // namespace tut

autotest/gdrivers/wmts.py

@@ -1245,7 +1245,8 @@ def test_wmts_20():
     <ServiceMetadataURL xlink:href="/vsimem/wmts_20.xml"/>
 </Capabilities>""")
 
-    ds = gdal.Open('WMTS:/vsimem/wmts_20.xml')
+    ds = gdal.OpenEx('WMTS:/vsimem/wmts_20.xml',
+                     open_options = ["CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX_LIMITS=YES"])
     assert ds is not None
     assert ds.RasterXSize == 512
     assert ds.RasterYSize == 256

doc/source/drivers/raster/wmts.rst

@@ -65,8 +65,7 @@ The WMTS driver can open :
 
       gdalinfo "WMTS:http://maps.wien.gv.at/wmts/1.0.0/WMTSCapabilities.xml,layer=lb"
 
--  the *WMTS:* prefix with open options URL (required), LAYER,
-   TILEMATRIXSET, TILEMATRIX, ZOOM_LEVEL, STYLE, EXTENDBEYONDDATELINE.
+-  the *WMTS:* prefix with open options
 
    ::
 
@@ -78,6 +77,65 @@ layer has more than one style or a tile matrix set, a list of
 subdatasets will be returned. If there is only one layer, it will be
 opened on the default style and the first tile matrix set listed.
 
+Open options
+------------
+
+The following open options are available:
+
+- **URL**: URL (or filename for local files) to GetCapabilities response document.
+  Required if not specified in the connection string (e.g if using "WMTS:" only)
+
+- **LAYER**: Layer identifier
+
+- **TILEMATRIXSET**: Tile Matrix Set identifier, which determines the CRS into
+  which the layer will be exposed. Must be one of the listed tile matrix
+  for the layer.
+
+- **TILEMATRIX**: Tile Matrix identifier. Must be one of the listed tile matrix of
+  the select tile matrix set for the layer. Mutually exclusive with ZOOM_LEVEL.
+  If not specified the last tile matrix, i.e. the one with the best resolution,
+  is selected.
+
+- **ZOOM_LEVEL**: Index of the maximum zoom level tile matrix to use for the
+  full resolution GDAL dataset (lower zoom levels will be used for overviews).
+  The first one (ie the one of lower resolution) is indexed 0.
+  Mutually exclusive with TILEMATRIX.
+  If not specified the last tile matrix, i.e. the one with the best resolution,
+  is selected.
+
+- **STYLE**: Style identifier. Must be one of layer.
+
+- **EXTENDBEYONDDATELINE** = YES/NO.  Whether to make the extent go over dateline
+  and warp tile requests. See ExtendBeyondDateLine parameter of the local service
+  description XML file described below for more details.
+
+- **EXTENT_METHOD** = AUTO/LAYER_BBOX/TILE_MATRIX_SET/MOST_PRECISE_TILE_MATRIX.
+  GDAL needs to retrieve an extent for the layer. Different sources are possible.
+  WGS84BoundingBox element at the Layer level, BoundingBox elements with potentially
+  several CRS at the Layer level, BoundingBox of the TileMatrixSet definitions
+  shared by all layers, and TileMatrixLimit definitions at the Layer level.
+  By default (AUTO), GDAL will try first with a WGS84BoundingBox/BoundingBox corresponding
+  to the CRS implied by the select TileMatrixSet. If not available, if will
+  fallback to a BoundingBox in another CRS and reproject it to the selected CRS.
+  If not available, it will fallback to the most precise tile matrix of the
+  selected TileMatrixSet and will clip it with the bounding box implied by the
+  most precise zoom level of the TileMatrixLimit of the layer.
+  If LAYER_BBOX is specified, only WGS84BoundingBox/BoundingBox elements are
+  considered.
+  If TILE_MATRIX_SET is specified, the BoundingBox element of the selected
+  TileMatrixSet will be used.
+  If MOST_PRECISE_TILE_MATRIX is specified, the implit extent of the
+  most precise tile matrix will be used.
+
+- **CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX** = YES/NO (GDAL >= 3.4.2).
+  Whether to use the implied bounds of the most precise TileMatrix to clip the
+  layer extent (defaults to NO if the layer bounding box is used, YES otherwise)
+
+- **CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX_LIMITS** = YES/NO (GDAL >= 3.4.2).
+  Whether to use the implied bounds of the most precise TileMatrixLimit to clip the
+  layer extent (defaults to NO if the layer bounding box is used, YES otherwise)
+
+
 Local service description XML file
 ----------------------------------
 

frmts/wmts/wmtsdataset.cpp

@@ -1503,6 +1503,8 @@ GDALDataset* WMTSDataset::Open(GDALOpenInfo* poOpenInfo)
         else if( EQUAL(pszExtentMethod, "MOST_PRECISE_TILE_MATRIX") )
             eExtentMethod = MOST_PRECISE_TILE_MATRIX;
 
+        bool bAOIFromLayer = false;
+
         // Use in priority layer bounding box expressed in the SRS of the TMS
         if( (!bHasAOI || bExtendBeyondDateLine) &&
             (eExtentMethod == AUTO || eExtentMethod == LAYER_BBOX) &&
@@ -1511,6 +1513,7 @@ GDALDataset* WMTSDataset::Open(GDALOpenInfo* poOpenInfo)
             if( !bHasAOI )
             {
                 sAOI = aoMapBoundingBox[oTMS.osSRS];
+                bAOIFromLayer = true;
                 bHasAOI = TRUE;
             }
 
@@ -1785,6 +1788,7 @@ GDALDataset* WMTSDataset::Open(GDALOpenInfo* poOpenInfo)
                                 sAOI.MaxY = std::max(std::max(dfY1, dfY2),
                                                     std::max(dfY3, dfY4));
                                 bHasAOI = TRUE;
+                                bAOIFromLayer = true;
                             }
                             delete poCT;
                         }
@@ -1798,7 +1802,7 @@ GDALDataset* WMTSDataset::Open(GDALOpenInfo* poOpenInfo)
         if( !bHasAOI && oTMS.bBoundingBoxValid &&
             (eExtentMethod == AUTO || eExtentMethod == TILE_MATRIX_SET) )
         {
-            CPLDebug("WMTS", "Using TMS bounding box");
+            CPLDebug("WMTS", "Using TMS bounding box as layer extent");
             sAOI = oTMS.sBoundingBox;
             bHasAOI = TRUE;
         }
@@ -1808,7 +1812,8 @@ GDALDataset* WMTSDataset::Open(GDALOpenInfo* poOpenInfo)
             (eExtentMethod == AUTO || eExtentMethod == MOST_PRECISE_TILE_MATRIX) )
         {
             const WMTSTileMatrix& oTM = oTMS.aoTM.back();
-            CPLDebug("WMTS", "Using TM level %s bounding box", oTM.osIdentifier.c_str() );
+            CPLDebug("WMTS", "Using TM level %s bounding box as layer extent",
+                     oTM.osIdentifier.c_str() );
 
             sAOI.MinX = oTM.dfTLX;
             sAOI.MaxY = oTM.dfTLY;
@@ -1826,41 +1831,79 @@ GDALDataset* WMTSDataset::Open(GDALOpenInfo* poOpenInfo)
             return nullptr;
         }
 
+        if( CPLTestBool(CSLFetchNameValueDef(
+                poOpenInfo->papszOpenOptions,
+                "CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX",
+                bAOIFromLayer ? "NO" : "YES")) )
         {
             // Clip with implied BoundingBox of the most precise TM
             // Useful for http://tileserver.maptiler.com/wmts
             const WMTSTileMatrix& oTM = oTMS.aoTM.back();
+            OGREnvelope sAOINew(sAOI);
 
             // For https://data.linz.govt.nz/services;key=XXXXXXXX/wmts/1.0.0/set/69/WMTSCapabilities.xml
             // only clip in Y since there's a warp over dateline.
             // Update: it sems that the content of the server has changed since
             // initial coding. So do X clipping in default mode.
             if( !bExtendBeyondDateLine )
             {
-                sAOI.MinX = std::max(sAOI.MinX, oTM.dfTLX);
-                sAOI.MaxX = std::min(sAOI.MaxX,
+                sAOINew.MinX = std::max(sAOI.MinX, oTM.dfTLX);
+                sAOINew.MaxX = std::min(sAOI.MaxX,
                     oTM.dfTLX +
                     oTM.nMatrixWidth  * oTM.dfPixelSize * oTM.nTileWidth);
             }
-            sAOI.MaxY = std::min(sAOI.MaxY, oTM.dfTLY);
-            sAOI.MinY =
+            sAOINew.MaxY = std::min(sAOI.MaxY, oTM.dfTLY);
+            sAOINew.MinY =
                 std::max(sAOI.MinY,
                          oTM.dfTLY -
                          oTM.nMatrixHeight * oTM.dfPixelSize * oTM.nTileHeight);
+            if( sAOI != sAOINew )
+            {
+                CPLDebug("WMTS",
+                         "Layer extent has been restricted from "
+                         "(%f,%f,%f,%f) to (%f,%f,%f,%f) using the "
+                         "implied bounding box of the most precise tile matrix. "
+                         "You may disable this by specifying the "
+                         "CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX open option "
+                         "to NO.",
+                         sAOI.MinX, sAOI.MinY, sAOI.MaxX, sAOI.MaxY,
+                         sAOINew.MinX, sAOINew.MinY, sAOINew.MaxX, sAOINew.MaxY);
+            }
+            sAOI = sAOINew;
         }
 
         // Clip with limits of most precise TM when available
+        if( CPLTestBool(CSLFetchNameValueDef(
+                poOpenInfo->papszOpenOptions,
+                "CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX_LIMITS",
+                bAOIFromLayer ? "NO" : "YES")) )
         {
             const WMTSTileMatrix& oTM = oTMS.aoTM.back();
             if( aoMapTileMatrixLimits.find(oTM.osIdentifier) != aoMapTileMatrixLimits.end() )
             {
+                OGREnvelope sAOINew(sAOI);
+
                 const WMTSTileMatrixLimits& oTMLimits = aoMapTileMatrixLimits[oTM.osIdentifier];
                 double dfTileWidthUnits = oTM.dfPixelSize * oTM.nTileWidth;
                 double dfTileHeightUnits = oTM.dfPixelSize * oTM.nTileHeight;
-                sAOI.MinX = std::max(sAOI.MinX, oTM.dfTLX + oTMLimits.nMinTileCol * dfTileWidthUnits);
-                sAOI.MaxY = std::min(sAOI.MaxY, oTM.dfTLY - oTMLimits.nMinTileRow * dfTileHeightUnits);
-                sAOI.MaxX = std::min(sAOI.MaxX, oTM.dfTLX + (oTMLimits.nMaxTileCol + 1) * dfTileWidthUnits);
-                sAOI.MinY = std::max(sAOI.MinY, oTM.dfTLY - (oTMLimits.nMaxTileRow + 1) * dfTileHeightUnits);
+                sAOINew.MinX = std::max(sAOI.MinX, oTM.dfTLX + oTMLimits.nMinTileCol * dfTileWidthUnits);
+                sAOINew.MaxY = std::min(sAOI.MaxY, oTM.dfTLY - oTMLimits.nMinTileRow * dfTileHeightUnits);
+                sAOINew.MaxX = std::min(sAOI.MaxX, oTM.dfTLX + (oTMLimits.nMaxTileCol + 1) * dfTileWidthUnits);
+                sAOINew.MinY = std::max(sAOI.MinY, oTM.dfTLY - (oTMLimits.nMaxTileRow + 1) * dfTileHeightUnits);
+
+                if( sAOI != sAOINew )
+                {
+                    CPLDebug("WMTS",
+                             "Layer extent has been restricted from "
+                             "(%f,%f,%f,%f) to (%f,%f,%f,%f) using the "
+                             "implied bounding box of the most precise tile matrix. "
+                             "You may disable this by specifying the "
+                             "CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX_LIMITS open option "
+                             "to NO.",
+                             sAOI.MinX, sAOI.MinY, sAOI.MaxX, sAOI.MaxY,
+                             sAOINew.MinX, sAOINew.MinY, sAOINew.MaxX, sAOINew.MaxY);
+                }
+                sAOI = sAOINew;
             }
         }
 
@@ -2287,6 +2330,8 @@ void GDALRegister_WMTS()
 "       <Value>TILE_MATRIX_SET</Value>"
 "       <Value>MOST_PRECISE_TILE_MATRIX</Value>"
 "  </Option>"
+"  <Option name='CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX' type='boolean' description='Whether to use the implied bounds of the most precise tile matrix to clip the layer extent (defaults to NO if layer bounding box is used, YES otherwise)'/>"
+"  <Option name='CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX_LIMITS' type='boolean' description='Whether to use the implied bounds of the most precise tile matrix limits to clip the layer extent (defaults to NO if layer bounding box is used, YES otherwise)'/>"
 "</OpenOptionList>");
 
     poDriver->pfnOpen = WMTSDataset::Open;

ogr/ogr_core.h

@@ -151,6 +151,21 @@ class CPL_DLL OGREnvelope
         return MinX <= other.MinX && MinY <= other.MinY &&
                MaxX >= other.MaxX && MaxY >= other.MaxY;
     }
+
+    /** Return whether the current rectangle is equal to the other rectangle */
+    bool operator== (const OGREnvelope& other) const
+    {
+        return MinX == other.MinX &&
+               MinY == other.MinY &&
+               MaxX == other.MaxX &&
+               MaxY == other.MaxY;
+    }
+
+    /** Return whether the current rectangle is not equal to the other rectangle */
+    bool operator!= (const OGREnvelope& other) const
+    {
+        return !(*this == other);
+    }
 };
 
 } // extern "C++"