OpenCV 地理空间栅格数据实战：用 GDAL 加载卫星影像与 SRTM 数字高程模型并生成地形热力图

OpenCV 地理空间栅格数据实战用 GDAL 加载卫星影像与 SRTM 数字高程模型并生成地形热力图【免费下载链接】opencvOpen Source Computer Vision Library项目地址: https://gitcode.com/GitHub_Trending/opencv31/opencv本篇指南围绕 OpenCV 官方教程doc/tutorials/app/raster_io_gdal.markdown展开讲解如何在不单独依赖 GDAL API 的前提下仅通过cv::imread配合IMREAD_LOAD_GDAL标志读取 GIS 栅格数据卫星影像、SRTM DEM并基于图像与 DEM 的四个角坐标建立经纬度到像素的映射完成逐像素海拔查询、地形热力图渲染与海平面上升 10/50/100 米淹没区叠加。读完后你将掌握 OpenCV GDAL 解码器的构建开关、读取标志组合、DEM 位深校验方法以及一份可直接运行的完整示例程序samples/cpp/tutorial_code/imgcodecs/GDAL_IO/gdal-image.cpp的实现细节。背景为什么 OpenCV 需要 GDAL地理空间栅格数据Geospatial Raster是地理信息系统GIS和摄影测量领域的核心产品形态典型用途是表示影像和数字高程模型Digital Elevation Model, DEM。GDALGeospatial Data Abstraction Library是加载 GIS 影像的标准库支持 GeoTIFF、SRTM、DTED、ENVI、FITS 等数十种栅格格式。OpenCV 的imgcodecs模块内置了一个 GDAL 解码器使得上述格式可以直接被cv::imread打开无需用户程序单独链接 GDAL API。该教程设定的目标可以归纳为五点使用 OpenCVimread加载卫星影像使用 OpenCVimread加载 SRTM 数字高程模型已知图像与 DEM 各自的角坐标时将 DEM 的海拔数据关联到图像上求出每个像素的海拔展示一个易于实现的 terrain heat map地形热力图展示 DEM 数据与正射校正ortho-rectified影像结合的一个基础用法。示例以旧金山的一张 GeoTIFF 影像和一张 SRTM DEM 为输入输出两部分结果一张地形热力图以及一张标注若海湾水位上升 10、50、100 米哪些区域会被淹没的叠加图。编译前提OpenCV 的 GDAL 支持是可选组件OpenCV 对 GDAL 的支持默认关闭。从顶层构建配置 CMakeLists.txt 可以看到OCV_OPTION(WITH_GDAL Include GDAL Support OFF VISIBLE_IF NOT ANDROID AND NOT IOS AND NOT XROS AND NOT WINRT VERIFY HAVE_GDAL)即WITH_GDAL默认为OFF且该选项在 Android、iOS、watchOS、WinRT 平台上不可见。要启用 GDAL 支持需要系统先安装 GDAL 开发库然后在 CMake 配置阶段显式打开WITH_GDALON。配置成功后构建摘要会打印出检测到的 GDAL 库与版本见 CMakeLists.txt 中的 status 输出逻辑if(WITH_GDAL OR HAVE_GDAL) status( GDAL: HAVE_GDAL THEN ${GDAL_LIBRARY} (ver ${GDAL_VERSION}) ELSE NO) endif()如果imread使用了 GDAL 标志但库未编译进去读取将失败返回空Mat因此在正式使用任何 GDAL 相关代码之前应先确认构建摘要中GDAL一行不是NO。用 GDAL 标志读取栅格数据GDAL 解码器的接入点在imread_内部。当标志位包含IMREAD_LOAD_GDAL时OpenCV 直接选用GdalDecoder跳过按文件签名查找常规解码器的流程见 modules/imgcodecs/src/loadsave.cpp#ifdef HAVE_GDAL if(flags ! IMREAD_UNCHANGED (flags IMREAD_LOAD_GDAL) IMREAD_LOAD_GDAL ){ decoder GdalDecoder().newDecoder(); }else{ #endif decoder findDecoder( filename );相关读取标志定义在 modules/imgcodecs/include/opencv2/imgcodecs.hpp 的ImreadModes枚举中与本文最相关的三个是标志值含义IMREAD_COLOR1总是转换为 3 通道 BGR 图像IMREAD_ANYDEPTH2输入是 16/32 位时按原始深度返回否则转为 8 位IMREAD_LOAD_GDAL8强制使用 GDAL 驱动加载图像加载卫星影像对普通影像只需IMREAD_LOAD_GDAL | IMREAD_COLORcv::Mat image cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );加载 DEM必须保留原始数值加载数字高程模型时每个像素的数值本身就是物理量海拔单位米不能被缩放或截断。教程用了一个直观的反例说明影像中一个以 double 表示值为 1 的像素和一个以 unsigned char 表示值为 255 的像素看起来完全一样但对地形数据而言像素值就是海拔高度。因此 DEM 必须加IMREAD_ANYDEPTH让 OpenCV 保留原生深度// load the dem model cv::Mat dem cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );如果你事先知道 DEM 的模型类型可以用Mat::type()或Mat::depth()做断言校验。NASA 或 DOD 的规范文档列出了各类高程模型的输入类型——SRTM 和 DTED 都是 signed short。示例程序里就做了这道安全门gdal-image.cpp// for sanity sake, make sure GDAL Loads it as a signed short if( dem.type() ! CV_16SC1 ){ throw std::runtime_error(DEM image type must be CV_16SC1); }这条断言并非拍脑袋从 GDAL 解码器的类型映射源码modules/imgcodecs/src/grfmt_gdal.cpp 中的gdal2opencv可以确认GDAL 的GDT_Int16会被映射为 OpenCV 的CV_16SC因此 SRTM.hgt文件signed short 存储在IMREAD_ANYDEPTH下读出来正好是CV_16SC1。若读出来的 Mat 类型不符说明深度转换发生了海拔数值就不可信了。同一文件中的range_cast函数grfmt_gdal.cpp展示了不带ANYDEPTH时会发生什么比如 16 位数据按 8 位输出时会执行value/256.0的整除截断——对影像这可能是可接受的视觉近似但对海拔 1234 米这样的数值则是灾难性的信息丢失。这也是教程强调 DEM 必须用IMREAD_ANYDEPTH的底层原因。grfmt_gdal.cpp中还可以看到解码器的一些实现细节构造时调用GDALAllRegister()注册全部 GDAL 驱动grfmt_gdal.cppcheckSignature通过文件头识别 NITF 与 DTED 标记grfmt_gdal.cpp逐行用band-RasterIO(GF_Read, ...)读入double扫描线再按目标深度写回 OpenCV Mat。此外imgcodecs.hpp 的注释提示使用IMREAD_LOAD_GDAL时ImageDecoder的返回值为 0 或 1因为 OpenCV 的 GDAL 解码器暂不支持多页multi-page读取。仓库自带的单元测试 modules/imgcodecs/test/test_gdal.cpp 验证了同样的标志组合img imread(path, cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH | cv::IMREAD_ANYCOLOR);该测试用 7 波段 float32 的 ENVI 与 FITS 文件断言读取出的像素精确值如EXPECT_EQ(101.125, img.atVecfloat,7(0,0)[0])说明 GDAL 路径对多波段、浮点深度的支持是可靠的。完整示例从影像到地形热力图与淹没区下面按 gdal-image.cpp 的程序结构逐段拆解。程序用法为gdal-image image_name dem_model_name第一步定义图像与 DEM 的角坐标示例假设图像是正射校正过的并预先给出四个角的经纬度GDAL 本身可以自动探测这些信息这里为演示目的硬编码// define the corner points // Note that GDAL library can natively determine this cv::Point2d tl( -122.441017, 37.815664 ); cv::Point2d tr( -122.370919, 37.815311 ); cv::Point2d bl( -122.441533, 37.747167 ); cv::Point2d br( -122.3715, 37.746814 ); // determine dem corners cv::Point2d dem_bl( -122.0, 38); cv::Point2d dem_tr( -123.0, 37);注意 DEM 只用了两个对角点dem_bl/dem_trSRTM 一个.hgt瓦片恰好覆盖 1°×1° 经纬度范围知道任意两角即可推出其余两角。第二步建立像素 → 世界 → DEM 像素的两次映射核心是两个线性插值函数pixel2world与world2demgdal-image.cpp。pixel2world把影像像素坐标换算成经纬度先在左右两条边上沿 y 方向插值出两条边界线上的点再沿 x 方向在两点间插值本质是一个双线性bilinear映射前提是图像正射校正后投影近似为经纬度的线性函数。cv::Point2d pixel2world( const int x, const int y, const cv::Size size ){ double rx (double)x / size.width; double ry (double)y / size.height; cv::Point2d rightSide lerp(tr, br, ry); cv::Point2d leftSide lerp(tl, bl, ry); return lerp( leftSide, rightSide, rx ); }world2dem则把经纬度反向映射回 DEM 像素坐标注释明确假设 DEM 数据是正射校正的cv::Point2d world2dem( cv::Point2d const coordinate, const cv::Size dem_size ){ double demRatioX ((dem_tr.x - coordinate.x)/(dem_tr.x - dem_bl.x)); double demRatioY 1-((dem_tr.y - coordinate.y)/(dem_tr.y - dem_bl.y)); cv::Point2d output; output.x demRatioX * dem_size.width; output.y demRatioY * dem_size.height; return output; }第三步逐像素取海拔并生成热力图主循环对图像每个像素做像素坐标 → 经纬度 → DEM 坐标 → 读取short型海拔值若 DEM 坐标越界则记为minElevation -10表示数据覆盖不到按低于海平面处理for( int y0; yimage.rows; y ){ for( int x0; ximage.cols; x ){ cv::Point2d coordinate pixel2world( x, y, image.size() ); cv::Point2d dem_coordinate world2dem( coordinate, dem.size() ); double dz; if( dem_coordinate.x 0 dem_coordinate.y 0 dem_coordinate.x dem.cols dem_coordinate.y dem.rows ){ dz dem.atshort(dem_coordinate); }else{ dz minElevation; } output_dem_flood.atcv::Vec3b(y,x) image.atcv::Vec3b(y,x); cv::Vec3b actualColor get_dem_color(dz); output_dem.atcv::Vec3b(y,x) actualColor; ... }}热力图的配色由一个颜色—海拔梯度表驱动gdal-image.cpp格式为pairVec3b 颜色, double 海拔从低到高推送color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 188, 154, 46), -1)); color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 110, 220, 110), 0.25)); color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 150, 250, 230), 20)); color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 160, 220, 200), 75)); color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 220, 190, 170), 100)); color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 250, 180, 140), 200));get_dem_color在相邻两个色标之间做线性插值低于最低色标返回最低色高于最高色标返回最高色否则定位所在区间并按比例t插值出Vec3bgdal-image.cpp。源码注释也提示这套渐变表非常适合外置成配置文件目前为了演示硬编码在程序里。第四步海平面上升淹没区叠加在拷贝原图的输出缓冲上按海拔分档叠加红、绿、蓝三种颜色分量直观表示水位上升 10 / 50 / 100 米后的受影响区域gdal-image.cpp// show effect of a 10 meter increase in ocean levels if( dz 10 ){ add_color( output_dem_flood.atcv::Vec3b(y,x), 90, 0, 0 ); } // show effect of a 50 meter increase in ocean levels else if( dz 50 ){ add_color( output_dem_flood.atcv::Vec3b(y,x), 0, 90, 0 ); } // show effect of a 100 meter increase in ocean levels else if( dz 100 ){ add_color( output_dem_flood.atcv::Vec3b(y,x), 0, 0, 90 ); }其中add_color对每个通道做了防溢出保护加色后限制在 0254。最后两个结果分别写出为heat-map.jpg和flooded.jpg。运行结果第一张是地形热力图低海拔区域海湾、谷底呈棕褐色随海拔升高过渡到绿、浅蓝、米黄到橙褐色。第二张是叠加在原卫星影像上的淹没分析——海湾沿岸低地泛红10 米档向内依次是绿50 米档与蓝100 米档的嵌套区域可以直观看出哪些街区和河谷在不同水位情景下会入水。数据准备方面输入影像可使用教程配图 gdal_output.jpg 对应的旧金山 GeoTIFFSRTM DEMN37W123.hgt从美国 USGS 的 SRTM v2.1 公开数据源获取教程原文给出了下载链接此处不再重复外部网址。实操要点与注意事项用 gdalinfo 查角坐标获取栅格角坐标的简便方法是命令行工具gdalinfo。对于已正射校正并带有投影信息的影像也可以直接从 USGS EarthExplorer 之类的数据门户读取。示例输出SRTM 瓦片N37W123.hgt$ gdalinfo N37W123.hgt Driver: SRTMHGT/SRTMHGT File Format Files: N37W123.hgt Size is 3601, 3601 Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, ... more output ... Corner Coordinates: Upper Left (-123.0001389, 38.0001389) (123d 0 0.50W, 38d 0 0.50N) Lower Left (-123.0001389, 36.9998611) (123d 0 0.50W, 36d5959.50N) Upper Right (-121.9998611, 38.0001389) (121d5959.50W, 38d 0 0.50N) Lower Right (-121.9998611, 36.9998611) (121d5959.50W, 36d5959.50N) Center (-122.5000000, 37.5000000) (122d30 0.00W, 37d30 0.00N) ... more output ...经纬度坐标在正式项目中应尽量避免教程的Notes部分特别指出地理坐标系Lat/Lon本质是球面坐标系把它当作笛卡尔坐标做线性插值在严格意义上是不正确的。本示例之所以直接使用经纬度是为了可读性且对演示目的足够精确。更严谨的做法是采用 UTMUniversal Transverse Mercator通用横轴墨卡托这类投影坐标系。其他工程性提示多波段/特殊格式从源码结构看GDAL 解码器按第一个波段的数据类型与波段数确定输出 Mat 类型并对GCI_PaletteIndex/GCI_GrayIndex/RGB 波段做颜色解释映射见 grfmt_gdal.cpp测试用例覆盖了 7 波段 float32 的 ENVI 文件说明多波段场景可用但输出通道数受IMREAD_ANYCOLOR/IMREAD_COLOR等标志约束。多页文件ImageDecoder文档注释表明 GDAL 解码器当前不支持多页读取imgcodecs.hpp遇到多页 TIFF 类文件时只能取第一页需要全量读取时应评估是否直接用 GDAL API。数值可信度DEM 场景中务必保留IMREAD_ANYDEPTH并配合类型断言如CV_16SC1否则 16 位海拔值会被按 256 倍缩放截断淹没区判断将完全失真。投影信息示例程序头部的注释明确说明——如果输入影像没有随附投影信息需要自己提供角坐标或改用完整 GDAL 驱动读取地理变换参数本例中的角点值是预先在文件顶部定义的。小结通过IMREAD_LOAD_GDAL与IMREAD_ANYDEPTH两个标志的组合OpenCV 让用户在纯 OpenCV 代码路径内完成了GIS 栅格 → Mat的桥接而本教程真正的技术含量在于其后的坐标关联用四角坐标做双向线性插值把 DEM 的海拔数值逐像素对齐到卫星影像上再经过色标渐变与分档叠色产出热力图与多档位淹没区分析图。这套像素 → 世界坐标 → DEM 像素的映射范式可以直接迁移到其他需要影像与高程/属性网格配准的场景生产化时建议将 UTM 投影、外置色标配置和投影元数据自动读取作为三个首要改进点。【免费下载链接】opencvOpen Source Computer Vision Library项目地址: https://gitcode.com/GitHub_Trending/opencv31/opencv创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考