使用C++/libCurl/Jsoncpp读取arcgis wmts 服务（restful模式）-白红宇

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

阅读量：6941 次

发布时间：2019-06-27

本文共 11221 字，大约阅读时间需要 37 分钟。

前言：

　　最近工作需要将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文件,内容如下:

<url>../data/world.tif</url> //底图

</image>

<url>http://localhost:6080/arcgis/rest/services/cj/MapServer</url>//瓦片服务

</image>

</map>

修改 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

另外附上编译好的osgEarth2.5