SRTM全球30米NASA高程DEM数据下载 msl高程

转载

hochie 2024-04-03 10:54:46

文章标签 c++ #include 数据图层 文章分类 机器学习人工智能

前言

osgearth_elevation示例，展示了如何通过点击地球获取不同定义下的高程数据。包括：MSL高程、HAE高程、EGM96高程。点击按钮，可以移除高程图层。

MSL高程：是mean sea level的简称，一般指平均海平面，是指水位高度等于观测结果平均值的平静的理想海面。观测时间范围不同，有不同概念的平均海平面，如日平均海平面、年平均海平面和多年平均海平面等等。一些验潮站常用18.6年或19年里每小时的观测值求出平均值，作为该站的平均海平面。简单记为：海平面高程。
HAE高程：WGS84坐标系下，定义椭球体，表示距椭圆体的高度的距离。简单记为：椭球体高程。
EFM96高程：EGM9模型是美国 BAI 推出的一种适用于全球范围，并综合利用现有全球大量重力数据所计算出来的高精度大地水准面模型。简单记为：米国自己定义的高程，好像在中国不太准确。

执行命令：

osgearth_elevationd.exe earth_image\world.earth

earth文件，添加了高程数据。否则无法获取高程值。

<map name="Globe" type="geocentric" version = "2">
	<!--此为全球影像图-->
	<image name="GlobeImage" driver="gdal">
		<url>./globe/globel.tif</url>
	</image>
	
	<!--全球高程图-->
	<heightfield name="GlobeHeightfiled" driver="gdal">
		<url>./heightfield/30m.tif</url>
	</heightfield>
	<!--文件缓存-->
	<options>
		<cache type="filesystem">
		  <path>./FileCache</path>
		</cache>
	</options>
</map>

执行效果

双击地形某处，或获取到高程等信息，并展示再ui面板上，同时会添加一个坐标模型。Scene graph interssection 是与地形交点处的高程。Query resolution 是分辨率。

SRTM全球30米NASA高程DEM数据下载 msl高程_图层

代码分析

#include <osgGA/StateSetManipulator>
#include <osgGA/GUIEventHandler>
#include <osgViewer/Viewer>
#include <osgViewer/ViewerEventHandlers>
#include <osgUtil/LineSegmentIntersector>
#include <osgEarth/MapNode>
#include <osgEarth/TerrainEngineNode>
#include <osgEarth/StringUtils>
#include <osgEarth/Terrain>
#include <osgEarth/VerticalDatum>
#include <osgEarthUtil/EarthManipulator>
#include <osgEarthUtil/Controls>
#include <osgEarthUtil/LatLongFormatter>
#include <osgEarthUtil/ExampleResources>
#include <osgEarthAnnotation/ModelNode>
#include <iomanip>

using namespace osgEarth;
using namespace osgEarth::Util;
using namespace osgEarth::Util::Controls;
using namespace osgEarth::Annotation;

// 创建一些全局静态对象
static MapNode*       s_mapNode     = 0L;	// 地图根结点
static LabelControl*  s_posLabel    = 0L;	// 显示经纬度
static LabelControl*  s_vdaLabel    = 0L;	// 显示坐标系geodetic(wgs84)
static LabelControl*  s_mslLabel    = 0L;	// 圆球之上的高度
static LabelControl*  s_haeLabel    = 0L;	// 椭球之上的高度
static LabelControl*  s_egm96Label  = 0L;	// EGM96的高度
static LabelControl*  s_mapLabel    = 0L;	// 场景图交点的高度
static LabelControl*  s_resLabel    = 0L;	// 输出分辨率结果
static ModelNode*     s_marker      = 0L;	// 坐标模型节点


// An event handler that will print out the elevation at the clicked point
// 一个事件处理程序，它将打印出单击点处的高程
struct QueryElevationHandler : public osgGA::GUIEventHandler 
{
    QueryElevationHandler()
        : _mouseDown( false ),
          _terrain  ( s_mapNode->getTerrain() )
    {
        _map = s_mapNode->getMap();
        _path.push_back( s_mapNode->getTerrainEngine() );
        _envelope = _map->getElevationPool()->createEnvelope(_map->getSRS(), 20u);// 20级LOD
    }

	// 重写update方法
    void update( float x, float y, osgViewer::View* view )
    {
        bool yes = false;

        // look under the mouse:查看鼠标下方
        osg::Vec3d world;
        osgUtil::LineSegmentIntersector::Intersections hits;// 保存与场景图 相交点的 std::multiset 表。
        if ( view->computeIntersections(x, y, hits) )// 计算与(x,y)相交的点，并存入hits
        {
            world = hits.begin()->getWorldIntersectPoint();

            // convert to map coords: 将交点转换到地形坐标
            GeoPoint mapPoint;
            mapPoint.fromWorld( _terrain->getSRS(), world );

            // do an elevation query:执行高程查询
            double query_resolution  = 0.0;  // max.全局并没有被用到
            double actual_resolution = 0.0;
            float elevation          = 0.0f;

			// 获取单个高程以及采样数据的分辨率。
            std::pair<float, float> result = _envelope->getElevationAndResolution(
                mapPoint.x(), mapPoint.y());

            elevation = result.first;			// msl高程
            actual_resolution = result.second;	// 分辨率

            if ( elevation != NO_DATA_VALUE )
            {
                // convert to geodetic to get the HAE:计算HAE下的高度
                mapPoint.z() = elevation;
                GeoPoint mapPointGeodetic( s_mapNode->getMapSRS()->getGeodeticSRS(), mapPoint );

                static LatLongFormatter s_f;

                s_posLabel->setText( Stringify()
                    << std::fixed << std::setprecision(2) 
                    << s_f.format(mapPointGeodetic.y(), true)
                    << ", " 
                    << s_f.format(mapPointGeodetic.x(), false) );

                if (s_mapNode->getMapSRS()->isGeographic())// 判断是否为地心投影
                {
                    double metersPerDegree = s_mapNode->getMapSRS()->getEllipsoid()->getRadiusEquator() / 360.0;
                    actual_resolution *= metersPerDegree * cos(osg::DegreesToRadians(mapPoint.y()));
                }

                s_mslLabel->setText( Stringify() << elevation << " m" );
                s_haeLabel->setText( Stringify() << mapPointGeodetic.z() << " m" );
                s_resLabel->setText( Stringify() << actual_resolution << " m" );// 分辨率

				// 计算egm96下的高度
                double egm96z = mapPoint.z();

                VerticalDatum::transform(
                    mapPointGeodetic.getSRS()->getVerticalDatum(),
                    VerticalDatum::get("egm96"),
                    mapPointGeodetic.y(),
                    mapPointGeodetic.x(),
                    egm96z);
                
                s_egm96Label->setText(Stringify() << egm96z << " m");

                yes = true;
            }

            // now get a normal ISECT HAE point.ESECT规定的HAE高程
            GeoPoint isectPoint;
            isectPoint.fromWorld( _terrain->getSRS()->getGeodeticSRS(), world );
            s_mapLabel->setText( Stringify() << isectPoint.alt() << " m");

            // and move the marker.添加模型
            s_marker->setPosition(mapPoint);

            // normal test.法线测试
            osg::Quat q;
            q.makeRotate(osg::Vec3(0,0,1), hits.begin()->getLocalIntersectNormal());
            s_marker->setLocalRotation(q);
        }

		// 根据yes 决定是否要更新标签
        if (!yes)
        {
            s_posLabel->setText( "-" );
            s_mslLabel->setText( "-" );
            s_haeLabel->setText( "-" );
            s_resLabel->setText( "-" );
            s_egm96Label->setText("-");
        }
    }

	// 接收左键点击、双击事件
    bool handle( const osgGA::GUIEventAdapter& ea, osgGA::GUIActionAdapter& aa )
    {
        if (ea.getEventType() == ea.DOUBLECLICK &&
            ea.getButton() == ea.LEFT_MOUSE_BUTTON)
        {
            osgViewer::View* view = static_cast<osgViewer::View*>(aa.asView());
			// 获取到点击事件，则更新坐标
            update( ea.getX(), ea.getY(), view );
            return true;
        }

        return false;
    }

    const Map*       _map;		// 地图
    const Terrain*   _terrain;	// 地形
    bool             _mouseDown;// 鼠标是否按下
    osg::NodePath    _path;		// 路径
    osg::ref_ptr<ElevationEnvelope> _envelope;
};

// 点击移除高程
struct ClickToRemoveElevation : public ControlEventHandler
{
    void onClick(Control*)
    {
        Map* map = s_mapNode->getMap();
        ElevationLayerVector layers;
        map->getLayers(layers);// 获取高程图层组
        map->beginUpdate();
        for (ElevationLayerVector::iterator i = layers.begin(); i != layers.end(); ++i) {
            map->removeLayer(i->get());// 移除高程数据
        }
        map->endUpdate();
    }
};


int main(int argc, char** argv)
{
    osg::ArgumentParser arguments(&argc,argv);

    osgViewer::Viewer viewer(arguments);

	// mapNode 节点
    s_mapNode = 0L;
    osg::Node* earthFile = MapNodeHelper().load(arguments, &viewer);
    if (earthFile)
        s_mapNode = MapNode::get(earthFile);

    if ( !s_mapNode )
    {
        OE_WARN << "Unable to load earth file." << std::endl;
        return -1;
    }

    osg::Group* root = new osg::Group();
    viewer.setSceneData( root );
    
    // install the programmable manipulator.
    viewer.setCameraManipulator( new osgEarth::Util::EarthManipulator() );

    // The MapNode will render the Map object in the scene graph.
    root->addChild( earthFile );

    // Make the readout: 设置布局
    Grid* grid = new Grid();
    grid->setBackColor(osg::Vec4(0,0,0,0.5));
    int r=0;
	// 第一行，都是标签，仅1个按钮——移除地形事件
    grid->setControl(0,r++,new LabelControl("Double-click to sample elevation", osg::Vec4(1,1,0,1)));
    grid->setControl(0,r++,new LabelControl("Coords (Lat, Long):"));
    grid->setControl(0,r++,new LabelControl("Map vertical datum:"));
    grid->setControl(0,r++,new LabelControl("Height above geoid:"));
    grid->setControl(0,r++,new LabelControl("Height above ellipsoid:"));
    grid->setControl(0,r++,new LabelControl("Scene graph intersection:"));
    grid->setControl(0,r++,new LabelControl("EGM96 elevation:"));
    grid->setControl(0,r++,new LabelControl("Query resolution:"));
    grid->setControl(0, r++, new ButtonControl("Click to remove all elevation data", new ClickToRemoveElevation()));

	// 第二行，标签
    r = 1;
    s_posLabel = grid->setControl(1,r++,new LabelControl(""));	// 显示经纬度
    s_vdaLabel = grid->setControl(1,r++,new LabelControl(""));	// 显示坐标系geodetic(wgs84)
    s_mslLabel = grid->setControl(1,r++,new LabelControl(""));	// 圆球之上的高度
    s_haeLabel = grid->setControl(1,r++,new LabelControl(""));	// 椭球之上的高度
    s_mapLabel = grid->setControl(1,r++,new LabelControl(""));	// 场景图交点的高度
    s_egm96Label = grid->setControl(1,r++,new LabelControl(""));// EGM96的高度
    s_resLabel = grid->setControl(1,r++,new LabelControl(""));	// 请求的结果

    
    Style markerStyle;
    markerStyle.getOrCreate<ModelSymbol>()->url()->setLiteral("D:/FreeXGIS/osgearth_gch/data/axes.osgt.64.scale");
    markerStyle.getOrCreate<ModelSymbol>()->autoScale() = true;
    s_marker = new ModelNode(s_mapNode, markerStyle);// 模型对象
    //s_marker->setMapNode( s_mapNode );
    //s_marker->setIconImage(osgDB::readImageFile("../data/placemark32.png"));
    s_marker->setDynamic(true);
    s_mapNode->addChild( s_marker );

	// 获取当前坐标系
    const SpatialReference* mapSRS = s_mapNode->getMapSRS();
    s_vdaLabel->setText( mapSRS->getVerticalDatum() ? 
        mapSRS->getVerticalDatum()->getName() : 
        Stringify() << "geodetic (" << mapSRS->getEllipsoid()->getName() << ")" );

	// 画布控件
    ControlCanvas* canvas = ControlCanvas::get(&viewer);
    canvas->addControl( grid );

    // An event handler that will respond to mouse clicks:
    viewer.addEventHandler( new QueryElevationHandler() );

    // add some stock OSG handlers:
    viewer.addEventHandler(new osgViewer::StatsHandler());
    viewer.addEventHandler(new osgViewer::WindowSizeHandler());
    viewer.addEventHandler(new osgGA::StateSetManipulator(viewer.getCamera()->getOrCreateStateSet()));

    return viewer.run();
}

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。