尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

OSG3D与osgEarth实现三维GIS水淹分析与长度量测

OSG3D与osgEarth实现三维GIS水淹分析与长度量测 简介OSG3D-master.zip是一套基于OpenSceneGraph与osgEarth的3D地理信息开发示例面向需要实现三维场景构建、空间测量与地形分析的OSG开发者尤其适合有一定OSG基础、希望快速上手osgEarth空间分析的用户。压缩包共88个文件以31个h头文件与29个cpp实现为主辅以earth场景配置、ui交互界面及vs工程构建文件整体仅66KB体量轻但结构完整便于快速定位核心代码。已有510人学习下载。示例明确分成OSG3D与OSGEarth两个工程分别演示点线面几何绘制以及长度、高度、面积测量和水淹分析功能其内部模块完整覆盖了洪水淹没模拟、两点与区域通视判断、Overlay叠加和全球网格测试等场景并结合数字高程模型直观展示实现思路。同时附带sln/vcxproj工程文件与说明文档为GIS二次开发、环境模拟和工程设计提供了可直接参考的代码框架。1. 从 OSG3D 到 osgEarth三维 GIS 场景里的水淹分析和长度量测水淹分析和长度测量在二维地图里做并不稀奇但你一旦把场景搬到三维地形上就会立刻碰到两类麻烦一是地形模型动辄上千万个三角形靠老办法逐片加载根本跑不动二是空间量测不能只看平面投影必须考虑高程起伏否则一段山路量出来和实际路径差一大截。OSG3DOpenSceneGraph和 osgEarth 的组合正好是解决这类问题的常用方案。OSG3D 提供高性能场景图负责几何节点、渲染状态和窗口交互osgEarth 则把全球影像和高程数据组织成动态调度的地形节点两者叠加后可以在同一个场景里完成水淹模拟和三维长度测量。这篇文章不是我写的项目实录而是基于这个标题最常见的技术路径做的一次完整梳理目的是让你看完后能直接用 C 在本地跑起来并理解每一步背后的几何原理。2. 用 OSG3D 几何节点搭起 osgEarth 地物框架加载地形与坐标投影2.1 osgEarth 地形初始化与 EPSG 坐标系在三维 GIS 场景里几何计算失败有八成是坐标系错误导致的。osgEarth 默认使用 EPSG:4326 经纬度坐标即 WGS84但渲染时它会内部投影到适合屏幕上显示的坐标系统。这里需要理解两个层面一是底图和高程数据以经纬度存储二是你在场景里添加几何体时必须使用 osgEarth 的 GeoPoint 进行坐标转换而不是直接把经纬度当成普通三维坐标塞进 osg::Geometry。坐标系用途常见场景EPSG:4326经纬度WGS84 基准原始定位数据、GeoPoint 参数EPSG:3857Web 墨卡托单位米导航地图、二维叠加自定义 UTM高斯投影单位米需要水平距离精确计算的局部区域初始化 osgEarth 场景的代码模板如下这个例子会在 114°E、40°N 附近加载一个影像图层并在该点放置一个黄色圆点作为地物标记#include osgViewer/Viewer #include osgEarth/MapNode #include osgEarth/Map #include osgEarth/EarthManipulator #include osgEarth/ImageLayer int main() { // 创建 osgEarth 地图并加载一个本地影像源 osg::ref_ptrosgEarth::Map map new osgEarth::Map(); osg::ref_ptrosgEarth::ImageLayer layer new osgEarth::ImageLayer( osgEarth::Config(image, osgEarth::Config(url, world.tif))); map-addLayer(layer.get()); osg::ref_ptrosgEarth::MapNode mapNode new osgEarth::MapNode(map.get()); // 用经纬度构造一个地物坐标 osgEarth::GeoPoint geoPoint( mapNode-getMapSRS()-getGeographicSRS(), 114.32, 40.04, 0.0); osg::Vec3d world; geoPoint.toWorld(world); // 创建 OSG3D 几何节点一个点 osg::ref_ptrosg::Geometry pointGeo new osg::Geometry(); osg::ref_ptrosg::Vec3Array vertices new osg::Vec3Array(); vertices-push_back(world); pointGeo-setVertexArray(vertices.get()); pointGeo-addPrimitiveSet(new osg::DrawArrays(osg::DrawArrays::POINTS, 0, 1)); osg::ref_ptrosg::Geode geode new osg::Geode(); geode-addDrawable(pointGeo.get()); mapNode-addChild(geode.get()); osgViewer::Viewer viewer; viewer.setSceneData(mapNode.get()); viewer.setCameraManipulator(new osgEarth::Util::EarthManipulator()); return viewer.run(); }这段代码的关键点在于GeoPoint::toWorld()它把经纬度转换为以场景原点为基础的三维世界坐标这样osg::Geometry里的顶点才能和地形渲染管线对齐。如果跳过这一步直接把(114.32, 40.04, 0.0)塞给Vec3Array你会看到一个点出现在渲染窗口中心附近但根本无法贴合地形。MapNode是所有 osgEarth 子节点的根相当于 OSG3D 场景图中的一个 Group但内部动态管理地形块和图层。你往这个 Group 下挂任何几何体都能和地形一起接受同样的裁剪和排序逻辑。2.2 用 osg::Geometry 创建可交互的几何标记osg::Geometry是 OSG3D 中最底层的几何绘制类它可以表示点、线、三角面也可组合多种图元。相比之下osgEarth 里的FeatureNode适合批量地理要素但它依赖样式解析动态修改变得笨重。做水淹区域和测量线段时我倾向于直接用osg::Geometry因为你可以完全控制顶点坐标、颜色数组和渲染顺序。一个要留意的地方是顶点顺序和环绕方向。沿地形表面构造多边形时顶点顺序必须是逆时针或顺时针统一否则背面剔除会导致多边形半透明或不可见。下面的代码演示如何动态生成一个三角面并绑定水深颜色osg::ref_ptrosg::Geometry surface new osg::Geometry(); // 四个顶点按逆时针排列 osg::ref_ptrosg::Vec3Array vts new osg::Vec3Array(); vts-push_back(osg::Vec3(0, 0, 10)); // 顶点 A vts-push_back(osg::Vec3(10, 0, 10)); // 顶点 B vts-push_back(osg::Vec3(10, 10, 10)); // 顶点 C vts-push_back(osg::Vec3(0, 10, 10)); // 顶点 D // 索引 osg::ref_ptrosg::DrawElementsUInt prim new osg::DrawElementsUInt( osg::DrawElementsUInt::TRIANGLE_FAN, 4); for (int i 0; i 4; i) prim-setElement(i, i); // 顶点颜色代表水深 0m、2m、4m、6m osg::ref_ptrosg::Vec4Array colors new osg::Vec4Array(); colors-push_back(osg::Vec4(0.2f, 0.4f, 1.0f, 0.6f)); colors-push_back(osg::Vec4(0.1f, 0.3f, 1.0f, 0.6f)); colors-push_back(osg::Vec4(0.0f, 0.2f, 1.0f, 0.6f)); colors-push_back(osg::Vec4(0.0f, 0.0f, 0.8f, 0.7f)); surface-setVertexArray(vts.get()); surface-setColorArray(colors.get(), osg::Array::BIND_PER_VERTEX); surface-addPrimitiveSet(prim.get()); surface-getOrCreateStateSet()-setMode(GL_BLEND, osg::StateAttribute::ON); surface-getOrCreateStateSet()-setRenderingHint(osg::StateSet::TRANSPARENT_BIN);BIND_PER_VERTEX意味着每个顶点都有自己的颜色插值会在三角形内部自动完成用来表现随水深渐变的效果非常合适。TRANSPARENT_BIN让这个几何体进入透明排序队列避免水面和地形互相遮挡产生闪烁。注意这里的坐标是直接写死的现实场景中应该由经纬度转换而来。3. 水淹模拟的几何实现动态水位面与地形淹没区间计算3.1 水位面生成与地形相交算法水淹模拟的核心问题可以描述为给定水位高度H找出地形表面所有高程低于H的点集合并把它们和地形融合成一个连续的视觉区域。实现思路有两种一种是用半透明颜色渲染不改变地形三角网另一种是直接构造一个与地形表面相交的多边形边界。后者更精准但需要三角形求交计算。osgEarth 提供osgEarth::ElevationQuery方法用来查询任意经纬度对应的地形高程。遍历目标区域内的网格点即可判定水下区域但直接遍历效率太低常见做法是把查询限制在矩形范围内然后按固定步长采样。下面是一个填充顶点数组的示例#include osgEarth/ElevationQuery osgEarth::ElevationQuery query(mapNode-getMap()); // 设定查询范围左下角 (114.30, 40.00)右上角 (114.36, 40.06) double x0 114.30, y0 40.00, x1 114.36, y1 40.06; double step 0.001; // 约 100 米一采样 int countX (x1 - x0) / step 1; int countY (y1 - y0) / step 1; osg::ref_ptrosg::Vec3Array waterVerts new osg::Vec3Array(); osg::ref_ptrconst osgEarth::SpatialReference srs mapNode-getMapSRS()-getGeographicSRS(); for (int i 0; i countX; i) { for (int j 0; j countY; j) { double lon x0 i * step; double lat y0 j * step; osgEarth::GeoPoint p(srs, lon, lat, 0.0); double elev 0.0; if (query.getElevation(p, elev, mapNode-getMapResourceCache())) { if (elev waterLevel) { // waterLevel 为当前水位高度 osg::Vec3d world; p.toWorld(world); waterVerts-push_back(world); } } } }getElevation是水淹判定中最常用的 API它接受GeoPoint和ResourceCache指针。返回值为true表示查询成功false表示该点在地图范围外。step参数直接影响结果精度和性能步长 0.001 度约等于 100 米用于宏观淹没范围足矣如果要做精细量算可以缩小到 0.00001 度但要保证MapNode的 LOD 足够高否则高程数据本身不会更新到网格层级。查询得到的顶点集只是一个无序点列直接绘制会形成碎片。你必须把它们转换为三角网网格。常用做法是构建 Delaunay 三角剖分或者利用地形本身的采样网格关系因为countX和countY是有序的可以按行列直接把相邻四个点组成两个三角形。这是最简单的结构化网格osg::ref_ptrosg::DrawElementsUInt triIndex new osg::DrawElementsUInt(GL_TRIANGLES); for (int j 0; j countY - 1; j) { for (int i 0; i countX - 1; i) { int idx j * countX i; triIndex-addElement(idx); triIndex-addElement(idx 1); triIndex-addElement(idx countX); triIndex-addElement(idx 1); triIndex-addElement(idx countX 1); triIndex-addElement(idx countX); } }这里假设所有点在网格中都有定义但淹没区域外的点并未加入waterVerts因此索引会错乱。我一般会预先创建一个std::vectorint标记每个网格点是否为水下点再通过映射表重新索引。这样避免了排序和哈希查找。3.2 淹没区域渲染与透明度参数得到三角网后你需要考虑视觉呈现方式。水淹区域一般用两种颜色混合表达浅水区带一点透明深水区全遮。可以沿用 2.2 里的BIND_PER_VERTEX思路在采样时计算每个点的水深depth waterLevel - elev并用归一化后的颜色插值输出。水深区间m颜色透明度0 - 1浅蓝 (0.3, 0.7, 1.0)0.551 - 3中蓝 (0.1, 0.4, 1.0)0.73 - 10深蓝 (0.0, 0.2, 0.8)0.8510 以上藏青 (0.0, 0.1, 0.4)1.0透明度值需要结合地形阴影调试。如果地形较暗半透明水面会变得难以辨识我习惯先给地形层关闭ambient或者调高天空光。另外osg::StateSet中需要关闭深度写入避免水面和地面产生深度竞争osg::StateSet* ss waterGeodetic-getOrCreateStateSet(); ss-setMode(GL_BLEND, osg::StateAttribute::ON); ss-setMode(GL_DEPTH_TEST, osg::StateAttribute::ON); ss-setAttributeAndModes(new osg::Depth(osg::Depth::LESS, 0.0, 1.0, false));osg::Depth构造函数最后一个参数false表示关闭深度写入这样水面不会影响后续地形渲染的深度缓冲。注意这里并不是标准动态水面反射而是静态水位高度下的一种模拟。如果要动画水位上升可以在每帧修改waterLevel并重新执行采样循环但这会带来一定性能开销后续第 5 章会说到动态更新的触发方式。4. 三维长度测量屏幕拾取、空间距离与贴合地形路径4.1 屏幕拾取与三维点获取长度测量用户点击地图上的两点引擎需要把屏幕像素坐标还原成地表三维坐标。最常用的工具是osgUtil::LineSegmentIntersector它从相机位置出发穿过鼠标位置发射线段求交。实现时需要在鼠标回调里取得窗口坐标并做以下转换#include osgUtil/LineSegmentIntersector #include osgViewer/Viewer #include osgGA/GUIEventHandler class PickAndMeasureHandler : public osgGA::GUIEventHandler { public: PickAndMeasureHandler(osgViewer::Viewer* viewer, osgEarth::MapNode* mapNode) : _viewer(viewer), _mapNode(mapNode) {} virtual bool handle(const osgGA::GUIEventAdapter ea, osgGA::GUIActionAdapter aa) { if (ea.getEventType() osgGA::GUIEventAdapter::RELEASE) { if (ea.getButton() 1) { // 左键 double x ea.getX(), y ea.getY(); osg::ref_ptrosgUtil::LineSegmentIntersector picker new osgUtil::LineSegmentIntersector( osgUtil::Intersector::PROJECTION, x, y); osgUtil::IntersectionVisitor iv(picker.get()); _viewer-getCamera()-accept(iv); if (picker-containsIntersections()) { const osgUtil::LineSegmentIntersector::Intersection result picker-getFirstIntersection(); osg::Vec3d world result.getWorldIntersectPoint(); accumulate(world); } } } return false; } private: osgViewer::Viewer* _viewer; osgEarth::MapNode* _mapNode; };LineSegmentIntersector在构造时默认是 WINDOW 坐标模式需要乘以视口宽高比例因此上面代码把ea.getX()和getY()直接作为投影坐标其实不严谨正确做法是调用ea.getXnormalized()和getYnormalized()或者先获得视口。为了减少兼容性问题我常用PROJECTION模式配合归一化坐标double nx ea.getXnormalized(); double ny ea.getYnormalized(); picker new osgUtil::LineSegmentIntersector(osgUtil::Intersector::PROJECTION, nx, ny);这样得到的world坐标是整个三维场景坐标它并不等于你直接读取的经纬度。要拿到地理坐标还需要做逆转换osgEarth::GeoPoint geo; geo.fromWorld(_mapNode-getMapSRS(), world); double lon geo.x(); double lat geo.y(); double elev geo.z();其中elev就是从当前地形表面拾取到的高程。注意这个高程不是原始 DEM 的高程因为地形图层可能叠加了多个数据源osgEarth 会综合计算结果。4.2 按地形高度修正测距结果两点之间的空间距离可以用欧几里得距离计算但这不符合 GIS 长度测量的习惯。测距目标通常是地表路径长度因此需要把相邻两个三维点投影到地形面上再沿连接线进行积分。最简单的做法是将路径细分成很多小段每一段都使用高程修正。如果你的测量范围不大可以先把world坐标转成平面投影坐标例如 UTM再用平面距离。代码如下double distanceInMeters(const osgEarth::GeoPoint p1, const osgEarth::GeoPoint p2, osgEarth::MapNode* mapNode) { // 将两个地形点转为同一投影坐标系 const osgEarth::SpatialReference* projSrs mapNode-getMapSRS()-getOrCreateProjectedSRS(); osgEarth::GeoPoint p1Proj p1.transform(projSrs); osgEarth::GeoPoint p2Proj p2.transform(projSrs); // 三维距离包含高程差 osg::Vec3d d(p2Proj.x() - p1Proj.x(), p2Proj.y() - p1Proj.y(), p2Proj.z() - p1Proj.z()); return d.length(); }投影坐标的单位是米所以length()返回的就是米。这里有个细节getOrCreateProjectedSRS()默认返回 UTM 带它会自动根据经纬度中心推断 UTM 分区。如果你的测量跨越多个 UTM 带结果会不稳定此时可改用等角割线投影或直接使用Geocentric坐标计算大地线弧长。下面是实际测量路径累积长度时的思路double total 0.0; for (size_t i 1; i measurePoints.size(); i) { total distanceInMeters(measurePoints[i-1], measurePoints[i], mapNode); }measurePoints是你通过多次点击收集到的GeoPoint数组。如果测量过程中点没有贴合地表而是停留在空中distanceInMeters的 z 分量过大导致结果失真。因此每新增一个点要强制让 y DEM 高程而不是 z 点击时的高程。可以用ElevationQuery重新采样double groundElev; query.getElevation(measurePoints.back(), groundElev, mapNode-getMapResourceCache()); measurePoints.back().z() groundElev;测量类型精度影响因素平面投影 高程忽略低地球曲率UTM 分带投影平面 高程直接相加中投影坐标系变形分段高程重采样 累加高采样密度地形 LOD在实际项目里我一般保留每个点击点的原始三维坐标同时保存对应的地形高程修正值这样可以随时切换显示“水平距离”和“实际距离”而不用重新点击一遍。5. 进阶动态水位更新时的 Geometry Dirty 标记与水面闪烁规避5.1 让水淹几何体跟随每帧更新如果你想让水位连续上升或下降就需要每帧更新水淹几何体的顶点位置。不少新手直接把Vec3Array里的顶点一个个赋值然后发现画面不刷新问题出在没有通知 OSG3D 几何体发生脏标记。osg::Geometry内部维护了显示列表顶点数组变了不调用相关 APIGPU 不会重新上传。正确做法void updateWaterVertices(osg::Geometry* waterGeom, const std::vectorosg::Vec3d vertices) { osg::Vec3Array* va static_castosg::Vec3Array*(waterGeom-getVertexArray()); for (size_t i 0; i vertices.size(); i) { (*va)[i] vertices[i]; } va-dirty(); waterGeom-dirtyDisplayList(); waterGeom-dirtyBound(); }dirty()告诉 GL 缓冲对象这个缓冲区内容已经改变dirtyDisplayList()让旧显示列表失效并重新生成dirtyBound()重新计算包围体避免视锥剔除提前剪掉水面。如果只改顶点不改颜色数组可以不改颜色。注意vertices.size()必须和原来的数组大小一致否则要重新设置顶点数组并标记VertexAttrib。5.2 水面与地形闪烁的处理闪烁多是深度冲突引起。水面和地形表面几乎共面时GPU 的浮点精度会让两边的像素随机交替。解决方法有三种第一种是给顶点高度增加一个微小偏移例如world.z() 0.3使水面浮在地形之上第二种是把地形层的深度纹理绑定到水面着色器里按深度比较结果决定是否丢弃片段第三种是使用 polygon offset在水面状态中启用osg::ref_ptrosg::PolygonOffset offset new osg::PolygonOffset(); offset-setFactor(-1.0f); offset-setUnits(-1.0f); ss-setAttributeAndModes(offset.get(), osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE);setFactor和setUnits负值是把多边形向外拉一点抵消共面冲突。实际调试时factor取 -1 到 -2units取 -1 到 -4 比较常见。5.3 用着色器绘制动态等高线验证水淹边界验证水淹范围是否准确的技巧之一是同时绘制一条水位线。以水位高度作为偏移把地形顶点 color 属性改为“是否淹没”你自己写一个简单的 GLSL 着色器用fract(worldPos.z * scale)实现等高线效果。下面是一个可运行的着色器片段// waterline.frag #version 120 uniform vec3 u_waterColor; varying vec3 v_worldPos; void main() { float height v_worldPos.z; float waterLevel 12.0; float width 0.2; if (abs(height - waterLevel) width) { gl_FragColor vec4(1.0, 0.8, 0.1, 1.0); } else if (height waterLevel) { gl_FragColor vec4(u_waterColor, 0.6); } else { discard; } }这个着色器思路很直接在片元阶段比较当前世界高度和水位阈值小于阈值的画水面等于阈值的画出一条亮线。顶点着色器必须将世界坐标传输出到v_worldPos否则无法判断。这种方案避免了 CPU 端重新生成几何体适合水位动态调整的调试环境。验证时你只需改变waterLevel这个 uniform就能观察一条黄线在坡度上移动直观检验采样步长是否合理。如果线呈锯齿状说明地形 LOD 不够或采样网格过粗而不是算法问题。本文还有配套的精品资源点击获取
返回列表