Mapbox地形切片及osgearth代码实现加载高程数据
本文参考：
osgearth加载mapbox在线高程数据
Mapbox Terrain-DEM地形切片原理浅析
本文分为三个部分：
1.mapbox部分地图数据简介
2.自建mapbox高程数据
3.使用osgearth加载mapbox高程数据
注：2.3. 部分都是别人博客内容，两个部分没有什么联系，但是使用步骤是前后关系。
一、MAPbox部分地图数据简介：
1.mapbox-terrain-dem-v1
官方介绍：
https://docs.mapbox.com/data/tilesets/reference/mapbox-terrain-dem-v1/
Mapbox Terrain-DEM v1 是 Mapbox 提供的光栅图块集，是一个全局海拔层。此图块集包含 PNG 图块的红色、绿色和蓝色通道中以米为单位的原始高度值，可以解码为以米为单位的原始高度。Mapbox Terrain-DEM 是 Mapbox Terrain-RGB v1 图块集的优化版本，具有一些更新的数据和一些压缩，以降低较低缩放级别的精度，从而制作更小、更快加载的图块。您可以将 Terrain-DEM 用于各种应用程序，包括视觉和分析，从样式地形斜坡和山影到为视频游戏生成 3D 地形网格。
2.mapbox-terrain-rgb-v1
官方介绍：
https://docs.mapbox.com/data/tilesets/reference/mapbox-terrain-rgb-v1/
Mapbox Terrain-RGB 是 Mapbox 提供的光栅瓦片集，包含以光栅 PNG 瓦片编码的全局高程数据，作为颜色值，可以解码为以米为单位的原始高度。您可以将 Terrain-RGB 用于各种视觉和分析应用，从设计地形斜坡和山影到为视频游戏生成 3D 地形网格。
================================================================
二、MAPbox数据切片及原理
1.Mapbox地形切片简介
Mapbox三维地形使用位移贴图将栅格卫星影像瓦片贴到相应的地形模型上，而mapbox的地形模型通过使用tif-RGB格式的全球数字高程数据切片为png或webp瓦片生成（tif-RGB格式数据比tif-dem轻量些），这些瓦片上每个像素的RGB值通过解码为以米为单位的原始高程值，然后使用mapbox渲染贴图形成三维地形。
那么如有生成离线地形需求时，就可以通过在NASA上按经纬度下载相应范围的geotiff格式的dem数据，然后将tif-dem（高程灰度图）转化/编码为tif-RGB格式，然后使用相应的切片工具进行切片（按照相应mapbox相应的规范组织瓦片），而这些处理在github上已经有许多开源的处理工具，这里推荐一个集成了上面所有处理的工具（根据dem数据生成地形切片的工具）,感兴趣的化也可以基于GDAL写一个。
对于Terrain-RGB的数据解码：Terrain-RGB 中每个颜色通道为 base-256，而三个颜色通道组成的为16,777,216 个值，这些值可以映射到 0.1 米的高度增量，从而实现三维地形应用所需的垂直精度，因此原来表示高程值是足够的。 mapbox收到图块后， 将需要获取各个像素的红色（R），绿色（G）和蓝色（B）值，然后根据下列公式解码为相应的高度值。
height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
2.Terrain-RGB与Terrain-DEM
对于mapbox在线地形服务Mapbox Terrain-DEM是Mapbox Terrain-RGB在线瓦片集的优化版本，主要为更新了高程数据与数据压缩，降低了低缩放级别下的精度等来实现体积更小的瓦片，而其它的比如：地区切片生成原理、数据解码等Mapbox Terrain-DEM与Mapbox Terrain-RGB都是相同的。
3.切片工具：
https://github.com/FreeGIS/dem2terrain
根据 DEM 数据生成地形切片工具，使用 NodeJS + GDAL（NodeBinding）开发制作。可用于用户自定义 DEM 高程数据源生产地形瓦片，以便局域网离线使用。
特点：
- 支持
mapbox和terrarium两种地形瓦片编码格式供mapboxgl使用，其中terrarium格式是tangram引擎的官方地形格式，tangram是另外一款开源的webgl二三维一体化的引擎； - 固定瓦片尺寸256，瓦片周围有1cell的buffer，即实际瓦片是258*258.
- 自动读取数据源的坐标系统，重编码输入的 DEM 栅格文件，并重投影至指定的坐标系4490、4326、3857，默认3857，然后生成瓦片；
- 支持适用于3857、4490、4326的地形切片生产；
- 内置了影像金字塔索引和多进程实现（暂未使用多线程），加速瓦片生成速度；
- 支持地形瓦片以文件目录或mbtiles两种格式存储；
- 命令行提供了瓦片生成的进图条提示，便于用户查看生成进度。
- 内置一些异常导致的临时文件清理工作。
三、OSGEARTH加载mapbox在线高程数据
mapbox的高程数据数据是rgba四通道的png格式栅格图像数据，而一般的高程数据HeightField是单通道的。
因此，在请求到源数据之后只需要做一次osg::Image -> osg::HeightField 的转换即可。
官网提供的转换公式：
height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
通过继承 TileSource 重写 createImage 和 createHeightField 即可实现功能。
头文件
#include <osgEarth/Common> #include <osgEarth/TileSource> #include <osgEarth/URI> using namespace osgEarth; class CgMapboxTerrainOptions : public TileSourceOptions { public: optional<URI>& url() { return _url; } const optional<URI>& url() const { return _url; } optional<bool>& invertY() { return _invertY; } const optional<bool>& invertY() const { return _invertY; } optional<std::string>& format() { return _format; } const optional<std::string>& format() const { return _format; } public: CgMapboxTerrainOptions(const TileSourceOptions& opt = TileSourceOptions()) : TileSourceOptions(opt) { setDriver("mapboxTer"); fromConfig(_conf); } CgMapboxTerrainOptions(const std::string& inUrl) : TileSourceOptions() { setDriver("mapboxTer"); fromConfig(_conf); url() = inUrl; } virtual ~CgMapboxTerrainOptions() { } public: Config getConfig() const { Config conf = TileSourceOptions::getConfig(); conf.updateIfSet("url", _url); conf.updateIfSet("format", _format); conf.updateIfSet("invert_y", _invertY); return conf; } protected: void mergeConfig(const Config& conf) { TileSourceOptions::mergeConfig(conf); fromConfig(conf); } private: void fromConfig(const Config& conf) { conf.getIfSet("url", _url); conf.getIfSet("format", _format); conf.getIfSet("invert_y", _invertY); } optional<URI> _url; optional<std::string> _format; optional<bool> _invertY; };
TileSource 的实现
通过继承 TileSource 重写 createImage 和 createHeightField，并自定义 TileSourceDriver 使用新的TileSource。
#include "CgMapboxTerrainOptions.h" #include <osgEarth/ImageUtils> #include <osgDB/FileNameUtils> #include <osgDB/FileUtils> #include <osgDB/Registry> #include <osgDB/ReadFile> #include <osgEarth/Registry> class MapBoxTerrainSource : public TileSource { public: MapBoxTerrainSource(const TileSourceOptions& options) :TileSource(options), _options(options), _rotate_iter(0u) { } Status initialize(const osgDB::Options* dbOptions) { _dbOptions = Registry::instance()->cloneOrCreateOptions(dbOptions); URI xyzURI = _options.url().value(); if (xyzURI.empty()) { return Status::Error("Fail: driver requires a valid \"url\" property"); } // driver requires a profile. if (!getProfile()) { return Status::Error("An explicit profile definition is required by the XYZ driver."); } _template = xyzURI.full(); _rotateStart = _template.find("["); _rotateEnd = _template.find("]"); if (_rotateStart != std::string::npos && _rotateEnd != std::string::npos && _rotateEnd - _rotateStart > 1) { _rotateString = _template.substr(_rotateStart, _rotateEnd - _rotateStart + 1); _rotateChoices = _template.substr(_rotateStart + 1, _rotateEnd - _rotateStart - 1); } _format = _options.format().isSet() ? *_options.format() : osgDB::getLowerCaseFileExtension(xyzURI.base()); return STATUS_OK; } osg::Image* createImage(const TileKey& key, ProgressCallback* progress) { unsigned x, y; key.getTileXY(x, y); if (_options.invertY() == true) { unsigned cols = 0, rows = 0; key.getProfile()->getNumTiles(key.getLevelOfDetail(), cols, rows); y = rows - y - 1; } std::string location = _template; // support OpenLayers template style: replaceIn(location, "${x}", Stringify() << x); replaceIn(location, "${y}", Stringify() << y); replaceIn(location, "${z}", Stringify() << key.getLevelOfDetail()); // failing that, legacy osgearth style: replaceIn(location, "{x}", Stringify() << x); replaceIn(location, "{y}", Stringify() << y); replaceIn(location, "{z}", Stringify() << key.getLevelOfDetail()); std::string cacheKey; if (!_rotateChoices.empty()) { cacheKey = location; unsigned index = (++_rotate_iter) % _rotateChoices.size(); replaceIn(location, _rotateString, Stringify() << _rotateChoices[index]); } URI uri(location, _options.url()->context()); if (!cacheKey.empty()) uri.setCacheKey(cacheKey); return uri.getImage(_dbOptions.get(), progress); } virtual std::string getExtension() const { return _format; } osg::HeightField* createHeightField(const TileKey& key, ProgressCallback* progress); private: const CgMapboxTerrainOptions _options; std::string _format; std::string _template; std::string _rotateChoices; std::string _rotateString; std::string::size_type _rotateStart, _rotateEnd; OpenThreads::Atomic _rotate_iter; osg::ref_ptr<osgDB::Options> _dbOptions; }; osg::HeightField* MapBoxTerrainSource::createHeightField(const TileKey& key, ProgressCallback* progress) { // height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1) // MapBox encoded elevation PNG. // https://www.mapbox.com/blog/terrain-rgb/ if (1/*_options.elevationEncoding().value() == "mapbox"*/) { if (getStatus().isError()) return 0L; osg::HeightField *hf = 0; osg::ref_ptr<osg::Image> image = createImage(key, progress); if (image.valid()) { // Allocate the heightfield. hf = new osg::HeightField(); hf->allocate(image->s(), image->t()); ImageUtils::PixelReader reader(image.get()); for (unsigned int c = 0; c < image->s(); c++) { for (unsigned int r = 0; r < image->t(); r++) { osg::Vec4 pixel = reader(c, r); pixel.r() *= 255.0; pixel.g() *= 255.0; pixel.b() *= 255.0; float h = -10000.0f + ((pixel.r() * 65536.0f + pixel.g() * 256.0f + pixel.b()) * 0.1f); hf->setHeight(c, r, h); } } } return hf; } else { return TileSource::createHeightField(key, progress); } } class MapBoxTerrainTileSourceDriver : public TileSourceDriver { public: MapBoxTerrainTileSourceDriver() { supportsExtension("osgearth_mapboxTer", "mapboxTer Driver"); } virtual const char* className() { return "mapboxTer Driver"; } virtual ReadResult readObject(const std::string& file_name, const Options* options) const { if (!acceptsExtension(osgDB::getLowerCaseFileExtension(file_name))) { return ReadResult::FILE_NOT_HANDLED; } return new MapBoxTerrainSource(getTileSourceOptions(options)); } }; REGISTER_OSGPLUGIN(osgearth_mapboxTer, MapBoxTerrainTileSourceDriver)
调用
CgMapboxTerrainOptions mbTerOpt; mbTerOpt.url() = "http://api.mapbox.com/v4/mapbox.terrain-rgb/{z}/{x}/{y}.pngraw?access_token=..."; mbTerOpt.profile()->namedProfile() = "spherical-mercator"; ElevationLayerOptions eoptions("mapboxEle", mbTerOpt); osg::ref_ptr<osgEarth::ElevationLayer> rpEleLayer = new osgEarth::ElevationLayer(eoptions); osg::ref_ptr<osgEarth::Map> rpMap = new osgEarth::Map(mapOpts); rpMap->addImageLayer(rpImagLayer.get()); rpMap->addElevationLayer(rpEleLayer.get());
重新定义配置选项 CgMapboxTerrainOptions 和 TileSourceDriver 后，可以将mapbox在线高程数据加载出来。
参考文章：
osgearth加载mapbox在线高程数据