地理空间栅格数据是地理信息系统和摄影测量中广泛使用的产品。栅格数据通常可以表示影像和数字高程模型（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++ 标准库
#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 );
 
// 确定DEM角点
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> output;
    for( int i=0; i<N; i++ ){
output[i] = (uchar)(((1-t)*minColor[i]) + (t * maxColor[i]));
    }
    return output;
}
 
/*
* 计算DEM颜色
*/
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 - color_range[x].second);
 
            break;
        }
    }
 
    // 颜色插值
    return lerp( color_range[idx].first, color_range[idx+1].first, t);
}
 
/*
* 给定像素坐标和输入图像的大小，计算像素位置
* 在DEM图像上。
*/
cv::Point2d world2dem( cv::Point2d const& coordinate, const cv::Size& dem_size ){
 
 
    // 将其与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将其加载为带符号短整型
    if( dem.type() != CV_16SC1 ){ throw std::runtime_error("DEM image type must be CV_16SC1"); }
 
    // 定义颜色范围以创建我们的输出DEM热力图
    // 配对格式（颜色，高程）；从低到高推入
    // 注意：这非常适合配置文件，但此处仅作为工作演示。
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;
 
    // 遍历图像中的每个像素，计算DEM点
    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，都是带符号的短整型。

注意事项

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

地理坐标系是球形坐标系，这意味着将其与笛卡尔数学结合使用在技术上是不正确的。本演示使用它们是为了提高可读性，并且其精度足以说明问题。更好的坐标系将是通用横轴墨卡托（Universal Transverse Mercator）。

查找角点坐标

查找图像角点坐标的一个简单方法是使用命令行工具 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 ...