前言：
最近工作需要将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">
<image name="t1" driver="gdal">
<url>../data/world.tif</url> //底图
</image>
<image name="t2" driver="arcgis">
<url>http://localhost:6080/arcgis/rest/services/cj/MapServer</url>//瓦片服务
</image>
</map>
修改 osgEarth中的application_osgearth_viewe工程，参数设置为
Debugging-->
Command Arguments : D:OSG_MAKEosgearth-2.5vs2010in est.earth -------test.earth文件位置
Woking Directory : ......in ------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工程中的 osgEarthsrcosgEarthDriversarcgisMapService.cpp
另外附上编译好的osgEarth2.5
https://pan.baidu.com/s/10F235iv7WHDZc24jhHMD7g 密码mwd9