OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.78k stars 2.5k forks source link

gdal_translate failed to generate the XML file of WMTS with token #5729

Closed Nature6T331 closed 2 years ago

Nature6T331 commented 2 years ago

The following URLs can be add in QGIS with XYZ Tiles and the map can be displayed(note ec899a50c7830ea2416ca182285236f3 is the token).

http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=ter&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TileMatrix={z}&TileRow={y}&TileCol={x}&tk=ec899a50c7830ea2416ca182285236f3

When I use gdal_ translate to generate wmts xml file, it failed.

gdal_translate "http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=ter&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TileMatrix={z}&TileRow={y}&TileCol={x}&tk=ec899a50c7830ea2416ca182285236f3" d:\tdtwmts.xml -of wmts

or gdal_translate "http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=ter&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TileMatrix={z}&TileRow={y}&TileCol={x}&tk=ec899a50c7830ea2416ca182285236f3" d:\tdtwmts.xml -of wmts --config gdal_http_tk ec899a50c7830ea2416ca182285236f3

Maybe GDAL cannot resolve parameter tk?

Nature6T331 commented 2 years ago

Here is the metadata URL: http://t0.tianditu.gov.cn/ter_w/wmts?request=GetCapabilities&service=wmts

Nature6T331 commented 2 years ago

May be we can define tk parameter in wmtsdataset.cpp?

if (CPLGetXMLValue(psXML, "tk", nullptr)) {
    CPLString optstr;
    optstr.Printf("tk=%s", CPLGetXMLValue(psXML, "tk", nullptr));
    http_request_opts = CSLAddString(http_request_opts, optstr.c_str());
}

WMTSAddOtherXML(psRoot, "tk", osOtherXML);

jratike80 commented 2 years ago

You should read https://gdal.org/drivers/raster/wmts.html again. You are using a GetTile template for creating the wmts xml file and that is your own invention, the documentation tell to use GetCapabilities. But unfortunately even if you do that part right you won't be able to read image data from that service because that service requires an extra &tk= parameter and the usage of tokens/api-keys is not standardized in the old generation of the OGC W*S services. For example the service you try to use http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetCapabilities&VERSION=1.0.0&tk=ec899a50c7830ea2416ca182285236f3 returns the GretTile url simply as <ows:Get xlink:href="http://t0.tianditu.com/ter_w/wmts?">

GDAL client follows the standard and drops the &tk= from GetTile requests.

For making the WMTS service to support standard clients it should modify the GetCapabilities so that the same token or api-key that was used in GetCapabilities would be added to all the links in the response. We do that with our services at the National Land Survey of Finland but that is certainly not common.

Maybe the GDAL WMTS xml schema could be enhanced to include a parameter "additional GetTile options" which would be added to each GetTile request.

Nature6T331 commented 2 years ago

Hi, jratike80, Thank you very much for your warmhearted answer. I probably knew what the problem was. The usage of tokens/api-keys is not standardized in WMTS driver, so &tk was dropped from requests. I have read https://gdal.org/drivers/raster/wmts.html many times. I use gdal_translate based on the example in the web page. Before that, my main purpose is to use GDAL open a WMTS map with token, which requires an XML file. I'm not familiar with the details of GDAL e especially the difference between GetTile and GetCapabilities.

Obviously, now I need to modify the source code of GDAL to achieve my goal. As you say, it should modify the GetCapabilities so that the same token or api-key that was used in GetCapabilities would be added to all the links in the response. Could you give me your service at the National Land Survey of Finland example?

I don't know how to enhance GDAL WMTS xml to include a parameter "additional GetTile options" which would be added to each GetTile request. Maybe rouault is too busy and has no time to improve it. Could you give me some suggestion?

In wmtsdataset.cpp, maybe we could add following code in function char** WMTSDataset::BuildHTTPRequestOpts(CPLString osOtherXML)

if (CPLGetXMLValue(psXML, "tk", nullptr)) { CPLString optstr; optstr.Printf("tk=%s", CPLGetXMLValue(psXML, "tk", nullptr)); http_request_opts = CSLAddString(http_request_opts, optstr.c_str()); } and WMTSAddOtherXML(psRoot, "tk", osOtherXML);

Where http_request_opts is used by CPLXMLNode* WMTSDataset::GetCapabilitiesResponse(const CPLString& osFilename, char** papszHTTPOptions)

More advanced is the ability to read generic additional HTTP parameter names, not only for tk. I have seen the source code of QGIS and can realize such a function, but I don't know how to realize it in GDAL.

jratike80 commented 2 years ago

Actually if you follow the documentation and use this syntax

gdal_translate "WMTS:http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetCapabilities&VERSION=1.0.0&tk=ec899a50c7830ea2416ca182285236f3,LAYER=ter" tdtwmts.xml -of wmts

then the generated tdtwmts.xml file will include the tk parameter

<GDAL_WMTS>
  <GetCapabilitiesUrl>http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&amp;REQUEST=GetCapabilities&amp;VERSION=1.0.0&amp;tk=ec899a50c7830ea2416ca182285236f3</GetCapabilitiesUrl>
  <Layer>ter</Layer>
  <Style>default</Style>
  <TileMatrixSet>w</TileMatrixSet>
  <DataWindow>
    <UpperLeftX>-20037512.19473402</UpperLeftX>
    <UpperLeftY>20037512.19473402</UpperLeftY>
    <LowerRightX>20037508.3427892</LowerRightX>
    <LowerRightY>-20037508.3427892</LowerRightY>
  </DataWindow>
  <BandsCount>4</BandsCount>
  <DataType>Byte</DataType>
  <Cache />
  <UnsafeSSL>true</UnsafeSSL>
  <ZeroBlockHttpCodes>204,404</ZeroBlockHttpCodes>
  <ZeroBlockOnServerException>true</ZeroBlockOnServerException>
</GDAL_WMTS>

With --debug on I can see errors like this but I do not know what they mean:

WMTS: Using tilematrix=14 (zoom level 13)
ERROR 1: Invalid dataset dimensions : 0 x 0

Anyway you have used gdal_translate in a wrong way and reported that as a bug. There may also be a bug but the recommended procedure is to ask help from the gdal-dev mailing list first.

Nature6T331 commented 2 years ago

Actually I have tried gdal_translate many ways. Once, I ran successfully and got a similar XML file. Now I use your syntex: gdal_translate "WMTS:http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetCapabilities&VERSION=1.0.0&tk=ec899a50c7830ea2416ca182285236f3,LAYER=ter" tdtwmts.xml -of wmts It returns back : ERROR 1: Invalid dataset dimensions : 0 x 0 Maybe sometimes it works.

I have asked the map server provider, they told me that the tk parameter was probably not passed. Judging from the source code, the tk parameter is not handled. How to ask help from the gdal-dev mailing list first?

jratike80 commented 2 years ago

When you created this issue you have seen this advice message:

Questions should go to the gdal-dev mailing list at https://lists.osgeo.org/mailman/listinfo/gdal-dev or other support forums. GitHub issues are for bug reports and suggestions for new features.

Go to https://lists.osgeo.org/mailman/listinfo/gdal-dev, create an account, and send mail.

Nature6T331 commented 2 years ago

Ok, I have send the email to gdal-dev@lists.osgeo.org. Thank you very much!

jratike80 commented 2 years ago

I did get a little bit further by using the suggested open option -oo CLIP_EXTENT_WITH_MOST_PRECISE_TILE_MATRIX=NO but now GDAL stops at error that the WMTS server is sending "Sorry, your request has been intercepted because it appears to be an attack."

The WGS84 bounding box in the GetCapabilities is wrong but I do not believe that it is the reason for the error

 <Layer>
            <ows:Title>ter</ows:Title>
            <ows:Abstract>ter</ows:Abstract>
            <ows:Identifier>ter</ows:Identifier>
            <ows:WGS84BoundingBox>
                <ows:LowerCorner>-20037508.3427892 -20037508.3427892</ows:LowerCorner>
                <ows:UpperCorner>20037508.3427892 20037508.3427892</ows:UpperCorner>
            </ows:WGS84BoundingBox>

This error may be the reason for the issue:

<TileMatrixSet>
            <ows:Identifier>w</ows:Identifier>
            <ows:SupportedCRS>urn:ogc:def:crs:EPSG::900913</ows:SupportedCRS>
            <TileMatrix>
                <ows:Identifier>1</ows:Identifier>
                <ScaleDenominator>2.958293554545656E8</ScaleDenominator>
                <TopLeftCorner>20037508.3427892 -20037508.3427892</TopLeftCorner>

First, it would be best to use EPSG:3857 instead of the old fake EPSG:900913 code but what bugs me is
`<TopLeftCorner>20037508.3427892 -20037508.3427892</TopLeftCorner>`.

The coordinates of the top left corner are in easting-northing order for EPSG:3857 and the first coordinate should be negative, the other positive. As the coordinates are now makes GDAL to create a request with these parameters <TileLevel>0</TileLevel> <TileX>-1</TileX> <TileY>-1</TileY> because it trusts the advertised tile matrix metadata. The WMTS server probably denies this request because it doen not have tiles with index (-1 -1).

I suggest to test GDAL and WMTS with a properly configured service first. The one used in the WMTS driver examples http://maps.wien.gv.at/wmts/1.0.0/WMTSCapabilities.xml should be good. That way you can learn how things should work, and after that you can try to find workarounds to deal with a WMTS service with erroneous metadata. The best option of course would be to get the service provider to fix their service.

Nature6T331 commented 2 years ago

I think the the service provider is OK. I have asked the map server provider, they told me that the tk parameter was probably not passed. It can open in QGIS with XYZ Tiles.

http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=ter&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TileMatrix={z}&TileRow={y}&TileCol={x}&tk=ec899a50c7830ea2416ca182285236f3

Also it can open in ArcMap use WMTS. In ArcMap ,we shoud customize a parameter named tk. By the way, if using QGIS open the xml file you genarated, it will show: appears to be an attack.

jratike80 commented 2 years ago

This line is definitely wrong even it may not be the reason for this issue with GDAL <TopLeftCorner>20037508.3427892 -20037508.3427892</TopLeftCorner> It may be that some clients do not interpret that line but fixing it should not make any harm to any client.

By the way, if using QGIS open the xml file you genarated, it will show: appears to be an attack. My theory is that GDAL interprets the wrong TopLeftCorner metadata and generates GetTile request by that information, that leads to negative tileX and tileY indexes, and server denies the request. The method for testing this theory is to fix the TopLeftCorner metadata that the server sends with GetCapabilities.

Nature6T331 commented 2 years ago

I have gone to https://lists.osgeo.org/mailman/listinfo/gdal-dev, created an account, and sent mail. There is no reply.

You can see the metadata: https://t0.tianditu.gov.cn/ter_w/wmts?request=GetCapabilities&VERSION=1.0.0&service=wmts

-20037508.3427892 -20037508.3427892 20037508.3427892 20037508.3427892 -20037508.3427892 -20037508.3427892 20037508.3427892 20037508.3427892 I have a test.xml like this: `` `https://t0.tianditu.gov.cn/ter_c/wmts?request=GetCapabilities&VERSION=1.0.0&service=wmts&tk=ec899a50c7830ea2416ca182285236f3` `ter` `` `c` `tiles` `` `-180.0` `90` ` 180` ` -90` `` `EPSG:4490` `4` `Byte` `` `true` `204,404` `true` `ec899a50c7830ea2416ca182285236f3` `` I used gdalinfo command gdalinfo C:\Users\Administrator\Desktop\test.xml It seems work well. It returns: Files: C:\Users\Administrator\Desktop\test.xml Size is 3963368, 1981684 Coordinate System is: GEOGCRS["China Geodetic Coordinate System 2000", DATUM["China 2000", ELLIPSOID["CGCS2000",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4490]] Data axis to CRS axis mapping: 2,1 Origin = (-180.000000000000000,90.000000000000000) Pixel Size = (0.000090831846703,-0.000090831846703) Metadata: ABSTRACT=ter TITLE=ter Image Structure Metadata: INTERLEAVE=PIXEL Corner Coordinates: Upper Left (-180.0000000, 90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N) Lower Left (-180.0000000, -90.0000173) (180d 0' 0.00"W, 90d 0' 0.06"S) Upper Right ( 180.0000346, 90.0000000) (180d 0' 0.12"E, 90d 0' 0.00"N) Lower Right ( 180.0000346, -90.0000173) (180d 0' 0.12"E, 90d 0' 0.06"S) Center ( 0.0000173, -0.0000087) ( 0d 0' 0.06"E, 0d 0' 0.03"S) Band 1 Block=128x128 Type=Byte, ColorInterp=Red Overviews: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242 Mask Flags: PER_DATASET ALPHA Overviews of mask band: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242 Band 2 Block=128x128 Type=Byte, ColorInterp=Green Overviews: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242 Mask Flags: PER_DATASET ALPHA Overviews of mask band: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242 Band 3 Block=128x128 Type=Byte, ColorInterp=Blue Overviews: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242 Mask Flags: PER_DATASET ALPHA Overviews of mask band: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242 Band 4 Block=128x128 Type=Byte, ColorInterp=Alpha Overviews: 1981684x990842, 990842x495421, 495421x247711, 247711x123855, 123855x61928, 61928x30964, 30964x15482, 15482x7741, 7741x3870, 3870x1935, 1935x968, 968x484, 484x242
jratike80 commented 2 years ago

I think I have not seen your mail coming to the list yet but perhaps it will appear soon. I know how to read the metadata and the error I am talking about is later in the document, here (and repeats for other TileMatrixes)

<TileMatrixSet>
            <ows:Identifier>w</ows:Identifier>
            <ows:SupportedCRS>urn:ogc:def:crs:EPSG::900913</ows:SupportedCRS>
            <TileMatrix>
                <ows:Identifier>1</ows:Identifier>
                <ScaleDenominator>2.958293554545656E8</ScaleDenominator>
                <TopLeftCorner>20037508.3427892 -20037508.3427892</TopLeftCorner>
                <TileWidth>256</TileWidth>
                <TileHeight>256</TileHeight>
                <MatrixWidth>2</MatrixWidth>
                <MatrixHeight>2</MatrixHeight>
            </TileMatrix>

It is true that gdalinfo works. Gdalinfo reads and interprets the GetCapabilities. What fails is when GDAL generates the GetTile request. It generates request based on the advertised TileIndexSet metadata. I tried to explain earlier how wrong metadata leads into requests with negative tileX and tileY values.

You can have a look at the definitions of the google3857 TileMatrixSet definitions on this server as a comparison http://maps.wien.gv.at/wmts/1.0.0/WMTSCapabilities.xml. There the coordinates of TopLeftCorner are it the right order <TopLeftCorner>-20037508.3428 20037508.3428</TopLeftCorner>.

rouault commented 2 years ago

I have gone to https://lists.osgeo.org/mailman/listinfo/gdal-dev, created an account, and sent mail.

check your spam folder for a confirmation email sent by the list

Nature6T331 commented 2 years ago

Hi rouault, There is no confirmation email in my spam folder. I have sent the email to gdal-dev@lists.osgeo.org again. Thank you!

Nature6T331 commented 2 years ago

Hi jratike80, It is true that gdalinfo works. So it means that there is no problem interpreting the parameter tk? Do you mean that there is something wrong in TileIndexSet metadata? But it can open in QGIS with XYZ Tiles. Also it can open in ArcMap use WMTS. In ArcMap ,we shoud customize a parameter named tk.

In the WMTS standard, TopLeftCorner is a character sequence describing the coordinates of the upper left corner of the TileMatrixSet, which is composed of coordinates x and y. In the geographical coordinate system, the order of longitude before latitude is not in line with international practice. Aviation and maritime departments usually expect latitude to be before longitude. In case of emergency, different coordinate display may lead to unsafe factors. Although no standard clearly stipulates that latitude must precede longitude, generally speaking, the order of latitude before longitude will be adopted.

TopLeftCorner order of common coordinate systems: 4326 YX 3857 XY 4490 YX EPSG:0 XY Plane coordinate system XY More infomation is here(You can use IE translate the Chinese) https://www.cnblogs.com/zhoulujun/p/15059582.html

jratike80 commented 2 years ago

The tile matrix set <ows:Identifier>w</ows:Identifier> that we talk about is using this CRS by GetCapabilities <ows:SupportedCRS>urn:ogc:def:crs:EPSG::900913</ows:SupportedCRS> EPSG:900913 does not exist but we know that it means EPSG:3857. Definition of EPSG:3857 can be checked from the EPSG registry https://epsg.org/crs/wkt/id/3857 ["Easting (X)",east],AXIS["Northing (Y)",north]

The axis order of the TopLeftCorner of TileMatrixSet must follow the corresponding CRS, so in this case it must be Easting-Northing. See the "Table 14 — Parts of TileMatrix data structure" of the WMTS standard (bolding added by me)

TopLeftCorner Position in CRS coordinates of the top-left corner of this tile matrix

Nature6T331 commented 2 years ago

EPSG:900913 can be defined as fake EPSG in GDAL DATA

Do you mean TopLeftCorner must be -20037508.3428 20037508.3428 rather then 20037508.3428 -20037508.3428?

Please see following URL to know two type of Tianditu Map's tile rules with different EPSG(4326,3857/900913) https://blog.csdn.net/mrbaolong/article/details/119766101

Google XYZ: Z represents zoom level, z = zoom; The origin of XY is in the upper left corner, X is from left to right, and Y is from top to bottom.

TMS: the standard of open source products. The definition of Z is the same as that of Google; The origin of XY is in the lower left corner, x from left to right and Y from bottom to top.

Quadtree: the coding specification used by Microsoft Bing Maps. The definition of Z is the same as that of Google. Tiles at the same level are not represented by two XY dimensions, but only an integer, which is subject to the quadtree coding rules.

Baidu(Chinese map) XYZ: Z starts from 1 and divides the map into four tiles at the highest level; The origin of XY is at the position of longitude 0 and latitude 0, x from left to right and Y from bottom to top.

Nature6T331 commented 2 years ago

EPSG:4490 https://t0.tianditu.gov.cn/ter_c/wmts?request=GetCapabilities&service=wmts TopLeftCorner is 90.0 -180.0 You mean it must be -90.0 180.0?

EPSG:900913 https://t0.tianditu.gov.cn/ter_w/wmts?request=GetCapabilities&service=wmts TopLeftCorner is 20037508.3428 -20037508.3428 You mean it must be -20037508.3428 20037508.3428?

Nature6T331 commented 2 years ago

@rouault I have CC @even.rouault@spatialys.com It seems I can't send email to gdal-dev@lists.osgeo.org?

jratike80 commented 2 years ago

EPSG:4490 https://epsg.org/crs_4490/China-Geodetic-Coordinate-System-2000.html is latitude longitude -> top-left corner is 90 -180 As mentioned earlier EPSG:3857 is Easting-Northing -> top-left corner is -20037508.3428 20037508.3428.

EPSG:900913 can be defined as fake EPSG in GDAL DATA

Yes, but because there is a real code EPSG:3857 it would be best to use that.

Because we deal now with a WMTS service the documentation about Google XYZ, TMS, Quadtree, and Baidu are irrelevant. In the WMTS standard the WMTS client must be able to get all the information that it needs for handling the tiling scheme from the GetCapabilities document. How to do it right both on the server side and on the client side is documented in the OGC WMTS 1.0 standard.

Because this certain service has problems with WMTS I would consider to try read the service with a GDAL WMS TMS minidriver instead https://gdal.org/drivers/raster/wms.html.

I must apologize that because of my other duties I can't make more trials to help you within a next week or so.

jratike80 commented 2 years ago

Sorry, but I had to test one more thing. I got it to almost work with the files in the attached zip. File local.xml is an edited copy of the GetCapabilities and gdallocal.xml is a WMTS definition file that is using this local copy. I tested this combination with this command: gdal_translate -of png -outsize 400 400 gdallocal.xml out.png --config CPL_CURL_VERBOSE YES I had the both files in the same directory where I fired the gdal_translate command. What happens is that GDAL creates GetTile requests that work fine when copied to a browser, here is one example: http://t0.tianditu.com/ter_w/wmts?tk=ec899a50c7830ea2416ca182285236f3&service=WMTS&request=GetTile&version=1.0.0&layer=ter&style=default&format=tiles&TileMatrixSet=w&TileMatrix=1&TileRow=0&TileCol=1 Request returns a tile that is a jpeg tile (format=tiles does not tell it) and looks like this

image

However, the gdal_translate command that I run from the command window does not work because the server returns HTTP 403 Forbidden. I do not know why the server blocks the request because the tk= token is in place. Perhaps the server needs some specific user client or cookies or something. But I believe that with this edited GetCapabilities file GDAL is requesting the right tiles.

gdallocal.zip

Nature6T331 commented 2 years ago

Thank you very much jratike80. I appreciate your help. Few people in the world can help me to solve this problem. This test means TopLeftCorner is not the problem. Also you say TopLeftCorner is 90.0 -180.0, so service for EPSG:4490 ( https://t0.tianditu.gov.cn/**ter_c**/wmts?request=GetCapabilities&service=wmts) is right. But this service also does not works by gdal_translate. If you can sure tk is in place by gdal_translate, so we must know why server returns HTTP 403 Forbidden.

Nature6T331 commented 2 years ago

I wish you can to try read the service with a GDAL WMS TMS minidriver instead. According to my cognitive ability, it may be difficult to make such an attempt. I'm not a professional programmer. My original purpose was to call this map. Thank you very much and look forward to your further testing

Nature6T331 commented 2 years ago

Map service providers insist that their services are OK because ArcGIS can be opened with WMTS. Maybe we can use fiddler to see why HTTP 403 Forbidden.

Nature6T331 commented 2 years ago

I have tried to use TMS with attached xml file in QGIS, the error message is testtms.zip :

IReadBlock failed at X offset 0, Y offset 0: GDALWMS: Unable to download block 0, 0. URL: HTTP status code: 401, error: (null). Add the HTTP status code to to ignore this error (see http://www.gdal.org/frmt_wms.html).,IReadBlock failed at X offset 0, Y offset 0: GDALWMS: Unable to download block 0, 0. URL: HTTP status code: 401, error: (null). Add the HTTP status code to to ignore this error (see http://www.gdal.org/frmt_wms.html).,IReadBlock failed at X offset 0, Y offset 0: GDALWMS: Unable to download block 0, 0. URL: HTTP status code: 401, error: (null). Add the HTTP status code to to ignore this error (see http://www.gdal.org/frmt_wms.html). IReadBlock failed at X offset 0, Y offset 0: GDALWMS: Unable to download block 0, 0. URL: HTTP status code: 401, error: (null). Add the HTTP status code to to ignore this error (see http://www.gdal.org/frmt_wms.html).

jratike80 commented 2 years ago

I found some time, I returned back to this challenge, and I tried to be logical.

Accept: text/html,application/xhtml+xml,application/xml;q=0.9,image/avif,image/webp,image/apng,*/*;q=0.8,application/signed-exchange;v=b3;q=0.9
Accept-Encoding: gzip, deflate
User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/101.0.4951.64 Safari/537.36

I used the local fixed GetCapabilities and gdallocal.xml file as I did in my previous test but now I enhanced the command to include the customized headers so that the whole command was

gdal_translate -of png -outsize 400 400 gdallocal.xml out.png --config curl_verbose on --config GDAL_HTTP_HEADER_FILE gdalheaders.txt
Input file size is 3963368, 3963368
0...10...20...30...40...50...60...70...80...90...100 - done.

The out.png was a valid png image showing a map of the whole World. I considered that my theory was correct.

Here is also a successful curl command curl -o out.png "http://t0.tianditu.com/ter_w/wmts?tk=ec899a50c7830ea2416ca182285236f3&service=WMTS&request=GetTile&version=1.0.0&layer=ter&style=default&format=tiles&TileMatrixSet=w&TileMatrix=1&TileRow=0&TileCol=1" -v -H "Accept: text/html,application/xhtml+xml,application/xml;q=0.9,image/avif,image/webp,image/apng,*/*;q=0.8,application/signed-exchange;v=b3;q=0.9" -H "Accept-Encoding: gzip, deflate" -H "User-Agent: Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/101.0.4951.64 Safari/537.36"

I had a try with AcrMAP with URL http://t0.tianditu.gov.cn/ter_w/wmts?request=GetCapabilities&service=wmts (&tk=added) and ArcMAP did connect the service, it showed the layer "ter" and I could add it into the map project but it did not draw a map for me. I did not try to debug ArcMAP and I have no plans to do that in the future.

You wrote that

This test means TopLeftCorner is not the problem.

but that is not true. In that test GDAL started to create valid http request with an edited TopLeftCorner. The failure in that test meant that the TopLeftCorner was not the only problem. I double checked that adding the headers is not enough by using the gdal_translate command above but this time with the original GetCapabilities with reversed EPSG:900913 TopLeftCorner x and y coordinates. The result was ERROR 1: Invalid dataset dimensions : 0 x 0.

The other service with "c" matrix set in https://t0.tianditu.gov.cn/**ter_c**/wmts?request=GetCapabilities&service=wmts should work if GDAL is told to use the custom headers but I did not test that.

Nature6T331 commented 2 years ago

@jratike80

I have some progress today.

First of all, I'm very sorry to tell you that the map service provider told me that there are two types of tokens: browser and server. The previous token can only be used in the browser. GDAL needs server token. I can send new token to you by email because it is inconvenient to make it public. Please tell me your email address.

Secondly, we tested the problem of tk parameter, and I also improved the code. It is necessary to add the following code to the file wmtsdataset.cpp:

        CPLString osToken = CPLURLGetValue(osCapabilitiesFilename, "tk");
        if (!osToken.empty())
            osURLTileTemplate = CPLURLAddKVP(osURLTileTemplate, "tk", osToken);

After recompilation, the XML file test can read the image. I can also send the XML file with the new token to your email. However, there are still some problems with the direct access by URL. Please test it.

Nature6T331 commented 2 years ago

Here is the file wmtsdataset.cpp. wmtsdataset.zip

Nature6T331 commented 2 years ago

Maybe we still should use tk=tk=ec899a50c7830ea2416ca182285236f3 The following URLs can be add in QGIS with XYZ Tiles and the map can be displayed(note ec899a50c7830ea2416ca182285236f3 is the token).

http://t0.tianditu.gov.cn/ter_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=ter&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TileMatrix={z}&TileRow={y}&TileCol={x}&tk=ec899a50c7830ea2416ca182285236f3

QGIS

jratike80 commented 2 years ago

I am sorry, I do have experience about GDAL as a user but I am not a programmer so I do not have any use for the wmtsdataset.cpp. But it seems that you have progressed well. When that service is used as a TMS layer in QGIS then QGIS is not really doing anything that is related to WMTS. The basic url is used as given and QGIS just inserts the values for {z}, {x}, and {y} . The WMTS capabilities are never read nor interpreted.

rouault commented 2 years ago

With the fix of https://github.com/OSGeo/gdal/pull/5821 , the following local WMTS XML file will work:

<GDAL_WMTS>
  <GetCapabilitiesUrl>https://t0.tianditu.gov.cn/ter_c/wmts?SERVICE=WMTS&amp;REQUEST=GetCapabilities&amp;VERSION=1.0.0</GetCapabilitiesUrl>
  <ExtraQueryParameters>tk=ec899a50c7830ea2416ca182285236f3</ExtraQueryParameters>
  <Layer>ter</Layer>
  <Style>default</Style>
  <TileMatrixSet>c</TileMatrixSet>
  <DataWindow>
    <UpperLeftX>-180</UpperLeftX>
    <UpperLeftY>90</UpperLeftY>
    <LowerRightX>180</LowerRightX>
    <LowerRightY>-90</LowerRightY>
  </DataWindow>
  <BandsCount>4</BandsCount>
  <DataType>Byte</DataType>
  <Cache />
  <UnsafeSSL>true</UnsafeSSL>
  <ZeroBlockHttpCodes>204,404</ZeroBlockHttpCodes>
  <ZeroBlockOnServerException>true</ZeroBlockOnServerException>
  <UserAgent>Mozilla/5.0</UserAgent>
  <Accept>*.*;q=0.8</Accept>
</GDAL_WMTS>

The key here to make things work is to specify ExtraQueryParameters, UserAgent and Accept with the above values (for the 2 later one, immitating what Firefox does, since the server seems to only accept browser clients). And also switch the X/Y order in DataWindow since the default one created by the driver doesn't work due to the server itself not correctly filling ows:WGS84BoundingBox and ows:BoundingBox (it should fill them with latitude, longitude order).

jratike80 commented 2 years ago

Well, OGC standards are not so straight always. The WGS84BoundingBox (mandatory in WMTS GetCapabilities) is defined in https://portal.ogc.org/files/?artifact_id=38867 "OGC Web Services Common Standard", in Table 34 — Parameters included in WGS84BoundingBox data type, and see it saying:

Ordered sequence of two double values in decimal degrees, with longitude before latitude

BoundingBox is defined in table 33 and it is

Ordered sequence of double values b

with a note

The number of axes included, and the order of these axes, shall be as specified by the referenced CRS.

The ExtraQueryParameters is an especially good enhancement, I have seen other services that require things like apikey, api-key, token.

rouault commented 2 years ago

you're right about WGS84BoundingBox that should be in long, lat order and BoundingBox follows the order of the CRS. I believe the driver behaves properly. So "https://t0.tianditu.gov.cn/ter_c/wmts?SERVICE=WMTS&REQUEST=GetCapabilities&VERSION=1.0.0&tk=ec899a50c7830ea2416ca182285236f3" output should be fixed to report

            <ows:BoundingBox>
                <ows:LowerCorner> -90.0 -180.0</ows:LowerCorner>
                <ows:UpperCorner>90.0 180.0</ows:UpperCorner>
            </ows:BoundingBox>

given that the implicit CRS is urn:ogc:def:crs:EPSG::4490

Nature6T331 commented 2 years ago

@rouault Thank you very much. This new code will be released in GDAL 3.5.1?