应用程序实用工具——使用 GDAL 读取地理空间光栅文件 OpenCV v4.8.0

上一个教程 ： 为我们的应用程序添加轨迹条

下一个教程 ： 使用 OpenCV 和相似性测量进行视频输入

原作者	Marvin Smith
兼容性	OpenCV >= 3.0

地理空间栅格数据是地理信息系统和摄影测量中大量使用的产品。栅格数据通常可以表示图像和数字高程模型（DEM）。加载 GIS 图像的标准库是地理数据抽象库（GDAL）。在本示例中，我们将展示使用 OpenCV 本地函数加载 GIS 栅格格式的技术。此外，我们还将举例说明 OpenCV 如何将这些数据用于新颖有趣的用途。

目标

本教程的主要目标是

• 如何使用 OpenCV imread 加载卫星图像。
• 如何使用 OpenCV imread - 加载 SRTM 数字高程模型
• 给定图像和 DEM 的角坐标，将高程数据与图像相关联，找出每个像素的高程。
• 展示一个简单易用的地形热图示例。
• 展示 DEM 数据与正射影像的基本用法。

为了实现这些目标，以下代码将数字高程模型和旧金山的 GeoTiff 图像作为输入。对图像和 DEM 数据进行处理后，将生成图像的地形热图，并标出当海湾水位上升 10 米、50 米和 100 米时会受影响的城市区域。

代码

/*
 * gdal_image.cpp -- 使用地理空间数据抽象库将 GIS 数据加载到 OpenCV 容器中
*/
// OpenCV 头文件
#include "opencv2/core.hpp
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp" */ // C++ 标准库
// C++ 标准库
#include <cmath
#include <iostream
#include <stdexcept
#include <vector
using namespace std；
// 定义角点
// 请注意，GDAL 库可以本地确定这一点
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 )；
// 确定演示角
cv::Point2d dem_bl( -122.0, 38)；
cv::Point2d dem_tr( -123.0, 37)；
// 热图颜色范围
std::vector<std::pair<cv::Vec3b,double> > color_range；
// 所有函数原型的列表
cv::Point2d lerp( const cv::Point2d&, const cv::Point2d&, const double& )；
cv::Vec3b get_dem_color( const double& )；
cv::Point2d world2dem( const cv::Point2d&, const cv::Size&)；
cv::Point2d pixel2world( const int&, const int&, const cv::Size& )；
void add_color( cv::Vec3b& pix, const uchar& b, const uchar& g, const uchar& r )；
/*
 * 线性插值
 * p1 - 点 1
 * p2 - 点 2
 * t - 从点 1 到点 2 的比率
*/
cv::Point2d lerp( cv::Point2d const& p1, cv::Point2d const& p2, const double& t ){
 return cv::Point2d( ((1-t)*p1.x) + (t*p2.x)、
 ((1-t)*p1.y) + (t*p2.y))；
}
/*
 * 插值颜色
*/
template <typename DATATYPE, int N>
cv::Vec<DATATYPE,N> lerp( cv::Vec<DATATYPE,N> const& minColor、
 cv::Vec<DATATYPE,N> const& maxColor、
 double const& t ){
 cv::Vec<DATATYPE,N> 输出；
 for( int i=0; i<N; i++ ){
 output[i] = (uchar)(((1-t)*minColor[i]) + (t * maxColor[i]))；
 }
 return output；
}
/*
 * 计算演示颜色
*/
cv::Vec3b get_dem_color( const double& elevation ){
 // 如果海拔高度低于最小值，返回最小值
 if( elevation < color_range[0].second ){
 return color_range[0].first；
 }
 // 如果海拔高度高于最大值，则返回最大值
 if( elevation > color_range.back().second ){
 return color_range.back().first；
 }
 // 否则，找到合适的起始索引
 int idx=0；
 double t = 0；
 for( int x=0; x<(int)(color_range.size()-1); x++ ){
 // 如果当前高度低于下一项，则使用当前的
 // 两种颜色作为我们的范围
 if( elevation < color_range[x+1].second ){
 idx=x；
 t = (color_range[x+1].second - elevation)/ (color_range[x+1].second - elevation)/ (color_range[x+1].second)
 (color_range[x+1].second - color_range[x].second)；
 break；
 }
 }
 // 插值颜色
 return lerp( color_range[idx].first, color_range[idx+1].first, t)；
}
/*
 * 给定像素坐标和输入图像的大小，计算像素在 DEM 图像上的位置。
 * 在 DEM 图像上的像素位置。
*/
cv::Point2d world2dem( cv::Point2d const& coordinate, const cv::Size& dem_size ){ // 将其与 DEM 点相关联。
 // 将其与 DEM 点相关联
 // 假设 dem 数据经过正交校正
 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；
}
/*
 * 将像素坐标转换为世界坐标
*/

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；
 // 计算每个坐标的 LERP
 cv::Point2d rightSide = lerp(tr, br, ry)；
 cv::Point2d leftSide = lerp(tl, bl, ry)；
 // 计算插值坐标的实际纬度/经度坐标
 return lerp( leftSide, rightSide, rx )；
}
/*
 * 为特定像素颜色值添加颜色
*/
void add_color( cv::Vec3b& pix, const uchar& b, const uchar& g, const uchar& r ){
 if( pix[0] + b < 255 && pix[0] + b >= 0 ）{ pix[0] += b; }
 if( pix[1] + g < 255 && pix[1] + g >= 0 ){ pix[1] += g; }
 if( pix[2] + r < 255 && pix[2] + r >= 0 ){ pix[2] += r; }
}
/*
 * 主函数
*/
int main( int argc, char* argv[] ){
 /*
 * 检查输入参数
 */
 if( argc < 3 ){
 cout << "usage： " << argv[0] << " <image_name> <dem_model_name>" << endl；
 return -1；
 }
 // 加载图像（注意，我们没有投影信息。您将
 // 需要自行加载或使用完整的 GDAL 驱动程序。这些值是预先定义的
 // 在此文件的顶部
 cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR )；
 // 加载 dem 模型
 cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH )；
 // 创建输出产品
 cv::Mat output_dem( image.size(), CV_8UC3 )；
 cv::Mat output_dem_flood( image.size(), CV_8UC3 )；
 // 为合理起见，确保 GDAL 将其加载为带符号的 short
 if( dem.type() != CV_16SC1 ){ throw std::runtime_error("DEM 图像类型必须是 CV_16SC1"); } }
 // 定义颜色范围以创建输出 DEM 热图
 // Pair format ( Color, elevation ); Push from low to high
 // 注意：这非常适合用于配置文件，但这里是用于工作演示。
 color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 188, 154, 46), -1))；
 color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 110, 220, 110), 0.25))；
 color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 150, 250, 230), 20)；
 color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 160, 220, 200), 75))；
 color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 220, 190, 170), 100))；
 color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 250, 180, 140), 200))；
 // 定义最低标高
 double minElevation = -10；
 // 遍历图像中的每个像素，计算演示点
 for( int y=0; y<image.rows; y++ ){
 for( int x=0; x<image.cols; x++ ){
 // 将像素坐标转换为纬/纬坐标
 cv::Point2d coordinate = pixel2world( x, y, image.size() )；
 // 根据纬度/纵向坐标计算 dem 图像像素坐标
 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.at<short>(dem_coordinate)；
 }else{
 dz = minElevation；
 }
 // 将像素值写入文件
 output_dem_flood.at<cv::Vec3b>(y,x) = image.at<cv::Vec3b>(y,x)；
 // 计算热图输出的颜色
 cv::Vec3b actualColor = get_dem_color(dz)；
 output_dem.at<cv::Vec3b>(y,x) = actualColor；
 // 显示海平面上升 10 米的影响
 if( dz < 10 ){
 add_color( output_dem_flood.at<cv::Vec3b>(y,x), 90, 0, 0 )；
 }
 // 显示海平面上升 50 米的影响
 else if( dz < 50 ){
 add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 90, 0 )；
 }
 // 显示海平面上升 100 米的影响
 else if( dz < 100 ){
 add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 0, 90 )；
 }
 }}
 // 打印热图
 cv::imwrite( "heat-map.jpg" , output_dem )；
 // 打印水浸效果图像
 cv::imwrite( "flooded.jpg", output_dem_flood)；
 return 0；
}

如何使用 GDAL 读取光栅数据

本演示使用默认的 OpenCV imread 函数。主要区别在于，为了强制 GDAL 加载图像，必须使用适当的标志。

 cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR )；

加载数字高程模型时，每个像素的实际数值至关重要，不能缩放或截断。例如，在图像数据中，一个像素表示为数值为 1 的二进制数，其外观与表示为数值为 255 的无符号字符的像素相同。对于地形数据，像素值代表海拔高度，单位为米。为了确保 OpenCV 保留原始值，请在 imread 中使用 GDAL 标志和 ANYDEPTH 标志。

 // 加载 dem 模型
 cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH )；

如果事先知道要加载的 DEM 模型的类型，那么使用断言或其他机制测试 Mat::type() 或 Mat::depth() 可能会比较安全。NASA 或 DOD 规范文档可以提供各种高程模型的输入类型。主要类型 SRTM 和 DTED 都是带符号的短线。

注释

通常应避免使用纬度/经度（地理）坐标

地理坐标系是一种球面坐标系，这意味着在笛卡尔数学中使用它们在技术上是不正确的。本演示使用这些坐标系是为了增加可读性，并且足够准确。更好的坐标系是通用横轴默卡托坐标系。

查找角坐标

查找图像角坐标的一个简单方法是使用命令行工具 gdalinfo。对于经过正射校正并包含投影信息的图像，可以使用 USGS EarthExplorer。

\f$> 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.50"W, 38d 0' 0.50"N)
 Lower Left (-123.0001389, 36.9998611) (123d 0' 0.50"W, 36d59'59.50"N)
 Upper Right (-121.9998611, 38.0001389) (121d59'59.50"W, 38d 0' 0.50"N)
 Lower Right (-121.9998611, 36.9998611) (121d59'59.50"W, 36d59'59.50"N)
 Center (-122.5000000, 37.5000000) (122d30' 0.00"W, 37d30' 0.00"N)
 ... more output ...

结果

以下是程序的输出结果。使用第一幅图像作为输入。有关 DEM 模型，请在此处下载位于 USGS 的 SRTM 文件。http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip

输入图像

热量图

热图叠加