用GDAL读取地理空间栅格文件
地理空间栅格数据是地理信息系统和摄影测量中使用最多的产品。光栅数据通常可以表示图像和数字高程模型(DEM)。加载GIS图像的标准库是地理数据抽象库(GDAL)。在本例中,我们将展示使用本机OpenCV函数加载GIS栅格格式的技术。另外,我们将会展示OpenCV如何将这些数据用于小而有趣的目的。
目标
本教程的主要目标:
- 如何使用OpenCV imread加载卫星图像。
- 如何使用OpenCV imread来加载SRTM数字高程模型
- 给定图像和DEM的角坐标,将高程数据校正到图像以找到每个像素的高程。
- 显示地形热图的基本,易于实现的示例。
- 显示基本使用的DEM数据加上正射校正图像。
为了实现这些目标,以下代码将数字高程模型以及旧金山的GeoTiff图像作为输入。处理图像和DEM数据,并生成图像的地形热图以及城市的水位上升10,50和100米时将受到影响的城市的标签区域。
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读取栅格数据
此演示使用默认的OpenCV imread功能。主要区别在于,为了强制GDAL加载映像,您必须使用相应的标志。
cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );
加载数字高程模型时,每个像素的实际数值是必不可少的,无法缩放或截断。例如,对于图像数据,值为1的表示为双倍的像素具有与表示为值为255的无符号字符的像素相同的外观。使用地形数据,像素值表示以米为单位的高程。为了确保OpenCV保留本机值,请在使用ANYDEPTH标志的imread中使用GDAL标志。
// 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都是签名短裤。
Notes
通常应避免纬度/经度(地理)坐标
地理坐标系是一个球面坐标系,这意味着使用它们与笛卡尔数学在技术上是不正确的。此演示使用它们来增加可读性,并且足够准确地指出了这一点。一个更好的协调系统将是通用横向墨卡托。
找到角坐标
找到图像的角坐标的一个简单方法是使用命令行工具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 ...
结果
下面是程序的输出。使用第一张图像作为输入。对于DEM模型,请下载位于USGS的SRTM文件。http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip
输入图像
热图
热图覆盖
更多建议: