用GDAL讀取地理空間柵格文件

2018-10-01 22:00 更新

地理空間柵格數(shù)據(jù)是地理信息系統(tǒng)和攝影測量中使用最多的產(chǎn)品。光柵數(shù)據(jù)通?？梢员硎緢D像和數(shù)字高程模型（DEM）。加載GIS圖像的標(biāo)準(zhǔn)庫是地理數(shù)據(jù)抽象庫（GDAL）。在本例中，我們將展示使用本機(jī)OpenCV函數(shù)加載GIS柵格格式的技術(shù)。另外，我們將會展示OpenCV如何將這些數(shù)據(jù)用于小而有趣的目的。

目標(biāo)

本教程的主要目標(biāo)：

  • 如何使用OpenCV imread加載衛(wèi)星圖像。
  • 如何使用OpenCV imread來加載SRTM數(shù)字高程模型
  • 給定圖像和DEM的角坐標(biāo)，將高程數(shù)據(jù)校正到圖像以找到每個像素的高程。
  • 顯示地形熱圖的基本，易于實(shí)現(xiàn)的示例。
  • 顯示基本使用的DEM數(shù)據(jù)加上正射校正圖像。

為了實(shí)現(xiàn)這些目標(biāo)，以下代碼將數(shù)字高程模型以及舊金山的GeoTiff圖像作為輸入。處理圖像和DEM數(shù)據(jù)，并生成圖像的地形熱圖以及城市的水位上升10,50和100米時將受到影響的城市的標(biāo)簽區(qū)域。

Code

/*
 * gdal_image.cpp -- Load GIS data into OpenCV Containers using the Geospatial Data Abstraction Library
*/
// OpenCV Headers
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
// C++ Standard Libraries
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
using namespace std;
// 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);
// range of the heat map colors
std::vector<std::pair<cv::Vec3b,double> > color_range;
// List of all function prototypes
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 );
/*
 * Linear Interpolation
 * p1 - Point 1
 * p2 - Point 2
 * t  - Ratio from Point 1 to Point 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));
}
/*
 * Interpolate Colors
*/
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;
}
/*
 * Compute the dem color
*/
cv::Vec3b get_dem_color( const double& elevation ){
    // if the elevation is below the minimum, return the minimum
    if( elevation < color_range[0].second ){
        return color_range[0].first;
    }
    // if the elevation is above the maximum, return the maximum
    if( elevation > color_range.back().second ){
        return color_range.back().first;
    }
    // otherwise, find the proper starting index
    int idx=0;
    double t = 0;
    for( int x=0; x<(int)(color_range.size()-1); x++ ){
        // if the current elevation is below the next item, then use the current
        // two colors as our range
        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;
        }
    }
    // interpolate the color
    return lerp( color_range[idx].first, color_range[idx+1].first, t);
}
/*
 * Given a pixel coordinate and the size of the input image, compute the pixel location
 * on the DEM image.
*/
cv::Point2d world2dem( cv::Point2d const& coordinate, const cv::Size& dem_size   ){
    // relate this to the dem points
    // ASSUMING THAT DEM DATA IS ORTHORECTIFIED
    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;
}
/*
 * Convert a pixel coordinate to world coordinates
*/
cv::Point2d pixel2world( const int& x, const int& y, const cv::Size& size ){
    // compute the ratio of the pixel location to its dimension
    double rx = (double)x / size.width;
    double ry = (double)y / size.height;
    // compute LERP of each coordinate
    cv::Point2d rightSide = lerp(tr, br, ry);
    cv::Point2d leftSide  = lerp(tl, bl, ry);
    // compute the actual Lat/Lon coordinate of the interpolated coordinate
    return lerp( leftSide, rightSide, rx );
}
/*
 * Add color to a specific pixel color value
*/
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; }
}
/*
 * Main Function
*/
int main( int argc, char* argv[] ){
    /*
     * Check input arguments
    */
    if( argc < 3 ){
        cout << "usage: " << argv[0] << " <image_name> <dem_model_name>" << endl;
        return -1;
    }
    // load the image (note that we don't have the projection information.  You will
    // need to load that yourself or use the full GDAL driver.  The values are pre-defined
    // at the top of this file
    cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );
    // load the dem model
    cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );
    // create our output products
    cv::Mat output_dem(   image.size(), CV_8UC3 );
    cv::Mat output_dem_flood(   image.size(), CV_8UC3 );
    // 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"); }
    // define the color range to create our output DEM heat map
    //  Pair format ( Color, elevation );  Push from low to high
    //  Note:  This would be perfect for a configuration file, but is here for a working demo.
    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));
    // define a minimum elevation
    double minElevation = -10;
    // iterate over each pixel in the image, computing the dem point
    for( int y=0; y<image.rows; y++ ){
    for( int x=0; x<image.cols; x++ ){
        // convert the pixel coordinate to lat/lon coordinates
        cv::Point2d coordinate = pixel2world( x, y, image.size() );
        // compute the dem image pixel coordinate from lat/lon
        cv::Point2d dem_coordinate = world2dem( coordinate, dem.size() );
        // extract the elevation
        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;
        }
        // write the pixel value to the file
        output_dem_flood.at<cv::Vec3b>(y,x) = image.at<cv::Vec3b>(y,x);
        // compute the color for the heat map output
        cv::Vec3b actualColor = get_dem_color(dz);
        output_dem.at<cv::Vec3b>(y,x) = actualColor;
        // show effect of a 10 meter increase in ocean levels
        if( dz < 10 ){
            add_color( output_dem_flood.at<cv::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.at<cv::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.at<cv::Vec3b>(y,x), 0, 0, 90 );
        }
    }}
    // print our heat map
    cv::imwrite( "heat-map.jpg"   ,  output_dem );
    // print the flooding effect image
    cv::imwrite( "flooded.jpg",  output_dem_flood);
    return 0;
}

如何使用GDAL讀取柵格數(shù)據(jù)

此演示使用默認(rèn)的OpenCV imread功能。主要區(qū)別在于，為了強(qiáng)制GDAL加載映像，您必須使用相應(yīng)的標(biāo)志。

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

加載數(shù)字高程模型時，每個像素的實(shí)際數(shù)值是必不可少的，無法縮放或截?cái)?。例如，對于圖像數(shù)據(jù)，值為1的表示為雙倍的像素具有與表示為值為255的無符號字符的像素相同的外觀。使用地形數(shù)據(jù)，像素值表示以米為單位的高程。為了確保OpenCV保留本機(jī)值，請?jiān)谑褂肁NYDEPTH標(biāo)志的imread中使用GDAL標(biāo)志。

    // load the dem model
    cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );

如果事先知道您正在加載的DEM模型的類型，那么使用斷言或其他機(jī)制來測試Mat :: type（）或Mat :: depth（）可能是一個安全的選擇。NASA或DOD規(guī)范文件可以為各種高程模型提供輸入類型。主要類型SRTM和DTED都是簽名短褲。

Notes

通常應(yīng)避免緯度/經(jīng)度（地理）坐標(biāo)

地理坐標(biāo)系是一個球面坐標(biāo)系，這意味著使用它們與笛卡爾數(shù)學(xué)在技術(shù)上是不正確的。此演示使用它們來增加可讀性，并且足夠準(zhǔn)確地指出了這一點(diǎn)。一個更好的協(xié)調(diào)系統(tǒng)將是通用橫向墨卡托。

找到角坐標(biāo)

找到圖像的角坐標(biāo)的一個簡單方法是使用命令行工具gdalinfo。對于正射校正并包含投影信息的圖像，可以使用USGS EarthExplorer。

\f{#content}gt; 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 ...

結(jié)果

下面是程序的輸出。使用第一張圖像作為輸入。對于DEM模型，請下載位于USGS的SRTM文件。http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip

GDAL

輸入圖像

GDAL

熱圖

GDAL

熱圖覆蓋

以上內(nèi)容是否對您有幫助：
在線筆記
App下載
App下載

掃描二維碼

下載編程獅App

公眾號
微信公眾號

編程獅公眾號