osgEarth的Rex引擎原理分析（六十一）如何提取出指定范围的高程网格

最新推荐文章于 2024-01-22 09:42:11 发布

hankern

最新推荐文章于 2024-01-22 09:42:11 发布

阅读量991

点赞数

分类专栏： osgEarth的Rex引擎原理分析 osgearth rex

本文链接：https://blog.csdn.net/hankern/article/details/104187543

版权

本文分析osgEarth的Rex引擎如何从DEM或IMG文件中提取指定范围的高程网格，涉及GDAL库的使用以及遥感影像的波段和波长概念。同时，列出了一系列关于osgEarth引擎内部工作原理的问题，涵盖从地球文件加载图层、地形瓦片管理到着色器和OpenGL资源的详细过程。

摘要由CSDN通过智能技术生成

目标：（五十五）中的问题132

如果从dem/img等文件中提取出指定范围的高程网格，需要用到gdal库进行读取。遥感影像一般包含至少一个波段，波段对应波长，比如红、绿、蓝等，也包含透明度。

osgEarthDrivers/gdal/ReaderWriterGDAL.cpp
osg::HeightField* createHeightField( const TileKey&        key,
                                         ProgressCallback*     progress)
    {
        if (key.getLevelOfDetail() > _maxDataLevel)
        {
            //OE_NOTICE << "Reached maximum data resolution key=" << key.getLevelOfDetail() << " max=" << _maxDataLevel <<  std::endl;
            return NULL;
        }

        GDAL_SCOPED_LOCK;

        int tileSize = getPixelsPerTile();

        //Allocate the heightfield
        osg::ref_ptr<osg::HeightField> hf = new osg::HeightField;
        hf->allocate(tileSize, tileSize);

        if (intersects(key))
        {
            //Get the meter extents of the tile
            double xmin, ymin, xmax, ymax;
            key.getExtent().getBounds(xmin, ymin, xmax, ymax);

            // Try to find a FLOAT band
            GDALRasterBand* band = findBandByDataType(_warpedDS, GDT_Float32);
            if (band == NULL)
            {
                // Just get first band
                band = _warpedDS->GetRasterBand(1);
            }

            if (_options.interpolation() == INTERP_NEAREST)
            {
                double colMin, colMax;
                double rowMin, rowMax;
                geoToPixel( xmin, ymin, colMin, rowMax );
                geoToPixel( xmax, ymax, colMax, rowMin );
                std::vector<float> buffer(tileSize * tileSize, NO_DATA_VALUE);

                int iColMin = floor(colMin);
                int iColMax = ceil(colMax);
                int iRowMin = floor(rowMin);
                int iRowMax = ceil(rowMax);
                int iNumCols = iColMax - iColMin + 1;
                int iNumRows = iRowMax - iRowMin + 1;

                int iWinColMin = max(0, iColMin);
                int iWinColMax = min(_warpedDS->GetRasterXSize()-1, iColMax);
                int iWinRowMin = max(0, iRowMin);
                int iWinRowMax = min(_warpedDS->GetRasterYSize()-1, iRowMax);
                int iNumWinCols = iWinColMax - iWinColMin + 1;
                int iNumWinRows = iWinRowMax - iWinRowMin + 1;

                int iBufColMin = osg::round((iWinColMin - iColMin) / double(iNumCols - 1) * (tileSize - 1));
                int iBufColMax = osg::round((iWinColMax - iColMin) / double(iNumCols - 1) * (tileSize - 1));
                int iBufRowMin = osg::round((iWinRowMin - iRowMin) / double(iNumRows - 1) * (tileSize - 1));
                int iBufRowMax = osg::round((iWinRowMax - iRowMin) / double(iNumRows - 1) * (tileSize - 1));
                int iNumBufCols = iBufColMax - iBufColMin + 1;
                int iNumBufRows = iBufRowMax - iBufRowMin + 1;

                int startOffset = iBufRowMin * tileSize + iBufColMin;
                int lineSpace = tileSize * sizeof(float);

                band->RasterIO(GF_Read, iWinColMin, iWinRowMin, iNumWinCols, iNumWinRows, &buffer[startOffset], iNumBufCols, iNumBufRows, GDT_Float32, 0, lineSpace);

                for (int r = 0, ir = tileSize - 1; r < tileSize; ++r, --ir)
                {
                    for (int c = 0; c < tileSize; ++c)
                    {
                        hf->setHeight(c, ir, buffer[r * tileSize + c]);
                    }
                }
            }
            else
            {
                double dx = (xmax