地理空間柵格資料是地理資訊系統和攝影測量中廣泛使用的產品。柵格資料通常可以表示影像和數字高程模型（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 ...