前言：

　　最近工作需要將arcgis的wmts服務接入我們的3DGis系統平臺，要求用戶只輸入一個rest模式的wmts服務地址，系統即可自動獲取並解析其元數據信息，生成wmts圖層，並渲染顯示。經過多種嘗試，最終通過參考修正osgEarth，獲得了我們需要的效果。過程中竟然花了3天編譯osgEarth，搞的很崩潰，還好最終搞定了。現將過程和收獲及教訓寫下！

正文：　　

　　開始計劃用libcurl獲取服務xml文檔，然後用libxml2進行解析，實際使用中發現 http://localhost:6080/arcgis/rest/services/cj/MapServer/WMTS/1.0.0/WMTSCapabilities.xml 使用 libcurl+libxml2解析http得到的字符串失敗，無法將字符串轉為xml文檔格式，自然也就沒法解析了，此路不通！！！

　　然後再另找方法，發現osgEarth有arcgis wmts 圖層，然後通過跟蹤調試osgEarh中的arcgis wmts代碼，初步能獲取 wmts元數據信息，但是osgEarth目前實現的arcgis wmts 服務只能是全球範圍，不能是局部一塊區域的，經過一番改造，終於達到了我們的目標，支持 地理坐標系 （全球及局部）

　　下載編譯osgEarth，測試其流程，參考 https://weibo.com/p/2304189447a8480102v2c2，編譯 osgEarth2.5 + osg3.0.1，根據我自己的wmts服務地址寫了一個test.earth文件,內容如下:

<map name="MyMap" type="geocentric" version="2">

修改 osgEarth中的application_osgearth_viewe工程，參數設置為

　　Debugging-->

　　　　Command Arguments : D:\OSG_MAKE\osgearth-2.5\vs2010\bin\test.earth　　　　　　-------test.earth文件位置

　　　　Woking Directory : ..\..\..\bin　　　　　　　　　　　　　　　　　　　　　　　　　　　------osgEarth運行目錄　

將osgdb_osgearth_arcgis和osgEarth工程的release版本的項目屬性做以下設置，以方便release模式下調試

　　　　C++ --->Optimization------>Disable(/Od)

　　　　linker ---> Debugging -----> Generate Debuf Info

運行調試（release版本），即可調試跟蹤到osgdb_osgearth_arcgis中的bool MapService::init( const URI& _uri, const osgDB::ReaderWriter::Options* options )函數中。

std::string sep = uri.full().find( "?" ) == std::string::npos ? "?" : "&";
std::string json_url = uri.full() + sep + std::string("f=pjson"); // request the data in JSON format

這兩句很關鍵！！！經過查詢我發現 arcgis可提供json格式的wmts元數據信息，url格式就是在rest服務地址後面加"?f=pjson"！！！！！

ReadResult r = URI(json_url).readString( options );
if ( r.failed() )
return setError( "Unable to read metadata from ArcGIS service" );

上面拯救readString跟蹤以後我們發現實際就是調用libcurl寫到std::stringstream中！！！調用過程如下：

　　　　　　　　osgEarth::TerrianLayer::initTileSource()
　　　　　　　　　　　　|
　　　　　　　　osgEarth::TileSource::startup(const osgDB::Options* options)
　　　　　　　　　　　　|
　　　　　　　　[osgdb_osgearth_argis.dll] ArcGISSource::initialize( const osgDB::Options* options)
　　　　　　　　　　　　 |
　　　　　　　　[osgdb_osgearth_argis.dll] MapService::init(const osgEarth::URI& _uri, const osgDB::Option* options)
　　　　　　{_baseURI="http://localhost:6080/arcgis/rest/services/cj/MapServer" _fullURI="http://localhost:6080/arcgis/rest/services/cj/MapServer" _cacheKey="" ...}
　　　　　　　　　　　　json_url : "http://localhost:6080/arcgis/rest/services/cj/MapServer?f=pjson"
　　　　　　　　　　　　|
　　　　　　　　ReadResult　　URI::readString(const osgDB::Options* dbOptions,ProgressCallback* progress ) const
　　　　　　　　　　　　 |
　　　　　　　　template<typename READ_FUNCTOR> ReadResult doRead(const URI& inputURI,const osgDB::Options* dbOptions,ProgressCallback* progress)
　　　　　　　　　　　　|
　　　　　　　　 result = reader.fromHTTP( uri.full(), remoteOptions.get(), progress );
　　　　　　　　　　　　|
　　　　　　　　ReadResult fromHTTP( const std::string& uri, const osgDB::Options* opt, ProgressCallback* p ) { return HTTPClient::readString(uri, opt, p); }
　　　　　　　　　　　　 |
　　　　　　　　ReadResult HTTPClient::readString(const std::string& location,const osgDB::Options* options,ProgressCallback* callback)
　　　　　　　　　　　　|
　　　　　　　　HTTPResponse response = this->doGet( location, options, callback );
　　　　　　　　　　　　|
　　　　　　　　HTTPResponse HTTPClient::doGet( const HTTPRequest& request, const osgDB::Options* options, ProgressCallback* callback) const !!!!!!核心代碼使用curl
得到的結果string如下:
"{ "currentVersion": 10.11, "serviceDescription": "", "mapName": "Layers", "description": "", ........

　// Read the profile. We are using "fullExtent"; perhaps an option to use "initialExtent" instead?
double xmin = 0.0;
double ymin = 0.0;
double xmax = 0.0;
double ymax = 0.0;
int srs = 0;

Json::Value fullExtentValue = doc["fullExtent"];
Json::Value extentValue = doc["extent"];
std::string srsValue;

// added a case for "extent" which can be removed if we want to fall back on initialExtent if fullExtent fails
if ( !fullExtentValue.empty() )
{
// if "fullExtent" exists .. use that
xmin = doc["fullExtent"].get("xmin", 0).asDouble();
ymin = doc["fullExtent"].get("ymin", 0).asDouble();
xmax = doc["fullExtent"].get("xmax", 0).asDouble();
ymax = doc["fullExtent"].get("ymax", 0).asDouble();
srs = doc["fullExtent"].get("spatialReference", Json::Value::null).get("wkid", 0).asInt();
}
else if( !extentValue.empty() )
{
// else if "extent" exists .. use that
xmin = doc["extent"].get("xmin", 0).asDouble();
ymin = doc["extent"].get("ymin", 0).asDouble();
xmax = doc["extent"].get("xmax", 0).asDouble();
ymax = doc["extent"].get("ymax", 0).asDouble();
srs = doc["extent"].get("spatialReference", Json::Value::null).get("wkid", 0).asInt();
}
else
{
// else "initialExtent" must exist ..
xmin = doc["initialExtent"].get("xmin", 0).asDouble();
ymin = doc["initialExtent"].get("ymin", 0).asDouble();
xmax = doc["initialExtent"].get("xmax", 0).asDouble();
ymax = doc["initialExtent"].get("ymax", 0).asDouble();
srs = doc["initialExtent"].get("spatialReference", Json::Value::null).get("wkid", 0).asInt();
}

//Assumes the SRS is going to be an EPSG code
std::stringstream sss;
sss << "epsg:" << srs;

std::string ssStr;
ssStr = sss.str();

if ( ! (xmax > xmin && ymax > ymin && srs != 0 ) )
{
GlbSetLastError(L"Map service does not define a full extent");
return false;
}

上面這段代碼是讀取wmts數據範圍和坐標系信息，下面在起始終止行列號的計算中會用到！

　 std::string format = "png";
int tile_rows = 256;
int tile_cols = 256;
int min_level = 25;
int max_level = 0;
int num_tiles_wide = 1;
int num_tiles_high = 1;

double origin_x = 0;
double origin_y = 0;

int start_x_col = 0;
int start_y_row = 0;
int end_x_col = 0;
int end_y_row = 0;

// Read the tiling schema
Json::Value j_tileinfo = doc["tileInfo"];
if ( !j_tileinfo.empty() )
{
// return setError( "Map service does not define a tiling schema" );

// TODO: what do we do if the width <> height?
tile_rows = j_tileinfo.get( "rows", 0 ).asInt();
tile_cols = j_tileinfo.get( "cols", 0 ).asInt();
if ( tile_rows <= 0 && tile_cols <= 0 )
{
GlbSetLastError(L"Map service tile size not specified");
return false;
}

format = j_tileinfo.get( "format", "" ).asString();
if ( format.empty() )
{
GlbSetLastError(L"Map service tile schema does not specify an image format");
return false;
}

Json::Value j_levels = j_tileinfo["lods"];
if ( j_levels.empty() )
{
GlbSetLastError(L"Map service tile schema contains no LODs");
return false;
}

Json::Value j_origin = j_tileinfo["origin"];
if ( j_origin.empty() )
{
GlbSetLastError(L"Map service tile schema contains no origin");
return false;
}

origin_x = j_origin.get( "x", 0).asDouble();
origin_y = j_origin.get( "y", 0).asDouble();

min_level = INT_MAX;
max_level = 0;
for( unsigned int i=0; i<j_levels.size(); i++ )
{
int level = j_levels[i].get( "level", -1 ).asInt();
if ( level >= 0 && level < min_level )
min_level = level;
if ( level >= 0 && level > max_level )
max_level = level;
}

if (j_levels.size() > 0)
{
int l = j_levels[0u].get("level", -1).asInt();
double res = j_levels[0u].get("resolution", 0.0).asDouble();
num_tiles_wide = (int)glb_round((xmax - xmin) / (res * tile_cols));
num_tiles_high = (int)glb_round((ymax - ymin) / (res * tile_rows));

//瓦片起始行列號(fixedTileLeftTopNumX、fixedTileLeftTopNumY)
//fixedTileLeftTopNumX = Math.floor((Math.abs(originX - minX))/(resolution*tileSize));
//fixedTileLeftTopNumY = Math.floor((Math.abs(originY - maxY))/(resolution*tileSize));
start_x_col = (int)floor(fabs(origin_x - xmin) / (res * tile_cols));
end_x_col = (int)floor(fabs(origin_x - xmax) / (res * tile_cols));
start_y_row = (int)floor(fabs(origin_y - ymin) / (res * tile_rows));
end_y_row = (int)floor(fabs(origin_y - ymax) / (res * tile_rows));
if (start_y_row > end_y_row)
swap(start_y_row,end_y_row);
if (start_x_col > end_x_col)
swap(start_x_col,end_x_col);

//In case the first level specified isn‘t level 0, compute the number of tiles at level 0
for (int i = 0; i < l; i++)
{
num_tiles_wide /= 2;
num_tiles_high /= 2;

start_x_col /= 2;
start_y_row /= 2;
end_x_col /= 2;
end_y_row /= 2;
}
//profile.setNumTilesWideAtLod0(num_tiles_wide);
//profile.setNumTilesHighAtLod0(num_tiles_high);

// 重新計算wmts瓦片真實坐標範圍
{// 地理坐標系
xmin = origin_x + start_x_col * res * tile_cols;
xmax = origin_x + (end_x_col+1) * res * tile_cols;
ymin = origin_y - (end_y_row+1) * res * tile_rows;
ymax = origin_y - start_y_row * res * tile_rows;
}
}
}

mpr_domMinLevel = min_level;
mpr_domMaxLevel = max_level;

mpr_domBlockSizeX = tile_cols;
mpr_domBlockSizeY = tile_rows;

mpr_pExtent = new CGlbExtent();
mpr_pExtent->Set(xmin,xmax,ymin,ymax);

mpr_startMinTileCol = start_x_col;
mpr_endMinTileCol = end_x_col;
mpr_startMinTileRow = start_y_row;
mpr_endMinTileRow = end_y_row;

上面這段是獲取tile分辨率，0級tile的行列數，0級tile的起始終止行列號以及wmts瓦片實際範圍。使用的是Jsoncpp庫作為json數據分析器，與c++兼容。

osgEarth中只計算了0級的行列數，我們實際應用中還需要用到 “起始和終止行列號” ，tile的實際範圍等信息（藍色部分）。

這裏有一個需要註意的地方 fullextent並不是wmts瓦片的實際範圍，而是比它們要小一些

所以需要重新wmts瓦片實際範圍，根據origin和res和tilewidth，tileheight信息可以計算出來。

然後根據此瓦片實際範圍計算出0級tile的真實坐標範圍，否則疊到地球上就會出現位置偏移，區域縮小問題！！

rest方式訪問tile圖片的規則如下：

http://192.168.1.211/arcgis/rest/services/cj/MapServer/level/row/col

wchar_t* buff = new wchar_t[MAX_PATH];
wsprintf(buff,L"%s/tile/%d/%d/%d",mpr_strUrl.c_str(),tileMatrix,tileRow,tileCol);

計算tile外包的方法

CGlbExtent* computeTileBound(glbInt32 tileMatrix, glbInt32 tileCol,glbInt32 tileRow)
{
// 天地圖的影像圖層最多到18級
if (tileMatrix > mpr_domMaxLevel) return NULL;

if (mpr_pExtent)
{
glbInt32 tileMatrixRows = abs(mpr_endMinTileRow - mpr_startMinTileRow + 1) * pow (2.0, tileMatrix-mpr_domMinLevel);
glbInt32 tileMatrixCols = abs(mpr_endMinTileCol - mpr_startMinTileCol + 1) * pow (2.0, tileMatrix-mpr_domMinLevel);
if (tileMatrixCols<=0 && tileMatrixRows<=0)
return NULL;
double tileMatrixColStep = mpr_pExtent->GetXWidth() / tileMatrixCols;
double tileMatrixRowStep = mpr_pExtent->GetYHeight() / tileMatrixRows;

int levelStartCol = mpr_startMinTileCol;
int levelStartRow = mpr_startMinTileRow;
for (int kk = mpr_domMinLevel+1; kk <= tileMatrix; kk++)
{
levelStartCol *= 2;
levelStartRow *= 2;
}

double minLongitude = mpr_pExtent->GetLeft() + (tileCol - levelStartCol) * tileMatrixColStep;
double maxLongitude = mpr_pExtent->GetLeft() + (tileCol - levelStartCol + 1) * tileMatrixColStep;
double minLatitude = mpr_pExtent->GetTop() - (tileRow - levelStartRow + 1) * tileMatrixRowStep;
double maxLatitude = mpr_pExtent->GetTop() - (tileRow - levelStartRow) * tileMatrixRowStep;

CGlbExtent* pe = new CGlbExtent();
pe->Set(minLongitude,maxLongitude,minLatitude,maxLatitude);
return pe;
}
}
return NULL;
}

經驗：

註意兩點: 我們測試使用的是arcgis 10.1版
1. arcgis wmts 服務名稱必須是 英文，不能是中文。否則讀不出json格式的元數據！！！！
2. arcgis wmts 服務發布順序為先分析、發布然後切片。否則會發現坐標系很可能從wgs84地理坐標系變成了wgs84投影坐標系！！

3. 使用http://localhost:6080/arcgis/rest/services/cj/MapServer?f=pjson方式獲取json格式的元數據。

4. 必須重新計算wmts瓦片的實際範圍。否則疊加到地球上後會發現會出現位置偏移和縮小問題！！！

參考文檔：

http://www.cnblogs.com/he-xiang/p/5679391.html ArcGisServer根據最大最小坐標換算瓦片行列號 中國小刀

https://blog.csdn.net/qq_25867649/article/details/52789467?locationNum=2 使用libcurl來下載文件

https://www.cnblogs.com/findumars/p/7252843.html C/C++使用libcurl庫發送http請求（get和post可以用於請求html信息，也可以請求xml和json等串）

osgEarth工程中的 osgEarth\src\osgEarthDrivers\arcgis\MapService.cpp

使用C++/libCurl/Jsoncpp讀取arcgis wmts 服務（restful模式）