深圳网站优化指导,西安专业做网站的公司,郴州网站开发,教人做美食视频网站返回:OpenCV系列文章目录#xff08;持续更新中......#xff09;
上一篇:OpenCV的周期性噪声去除滤波器(70)
下一篇 :OpenCV系列文章目录#xff08;持续更新中......#xff09; 目录
目标
代码#xff1a;
解释#xff1a;
如何使用 GDAL 读取栅格数据
注意
…返回:OpenCV系列文章目录持续更新中......
上一篇:OpenCV的周期性噪声去除滤波器(70)
下一篇 :OpenCV系列文章目录持续更新中...... 目录
目标
代码
解释
如何使用 GDAL 读取栅格数据
注意
通常应避免使用经度/纬度地理坐标
查找拐角坐标
结果
地理空间栅格数据是地理信息系统和摄影测量中大量使用的产品。栅格数据通常可以表示影像和数字高程模型 DEM。用于加载 GIS 影像的标准库是地理数据抽象库 GDAL。在此示例中我们将展示使用本机 OpenCV 函数加载 GIS 栅格格式的技术。此外我们将展示一些示例说明OpenCV如何将这些数据用于新颖有趣的目的。
目标
本教程的主要目标
如何使用 OpenCV imread 加载卫星图像。如何使用 OpenCV imread 加载 SRTM 数字高程模型给定图像和 DEM 的角坐标将高程数据与图像相关联以查找每个像素的高程。显示一个基本、易于实施的地形热图示例。显示 DEM 数据与正射校正影像的基本用法。
为了实现这些目标以下代码采用数字高程模型以及旧金山的 GeoTiff 图像作为输入。对影像和 DEM 数据进行处理并生成影像的地形热图并标注在海湾水位上升 10、50 和 100 米时将受到影响的城市区域。
代码如下
/** 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 vectorusing 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::vectorstd::paircv::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::VecDATATYPE,N lerp( cv::VecDATATYPE,N const minColor,cv::VecDATATYPE,N const maxColor,double const t ){cv::VecDATATYPE,N output;for( int i0; iN; 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 minimumif( elevation color_range[0].second ){return color_range[0].first;}// if the elevation is above the maximum, return the maximumif( elevation color_range.back().second ){return color_range.back().first;}// otherwise, find the proper starting indexint idx0;double t 0;for( int x0; x(int)(color_range.size()-1); x ){// if the current elevation is below the next item, then use the current// two colors as our rangeif( elevation color_range[x1].second ){idxx;t (color_range[x1].second - elevation)/(color_range[x1].second - color_range[x].second);break;}}// interpolate the colorreturn lerp( color_range[idx].first, color_range[idx1].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 ORTHORECTIFIEDdouble 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 dimensiondouble rx (double)x / size.width;double ry (double)y / size.height;// compute LERP of each coordinatecv::Point2d rightSide lerp(tr, br, ry);cv::Point2d leftSide lerp(tl, bl, ry);// compute the actual Lat/Lon coordinate of the interpolated coordinatereturn 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 dont 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 filecv::Mat image cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );// load the dem modelcv::Mat dem cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );// create our output productscv::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 shortif( 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::paircv::Vec3b,double(cv::Vec3b( 188, 154, 46), -1));color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 110, 220, 110), 0.25));color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 150, 250, 230), 20));color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 160, 220, 200), 75));color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 220, 190, 170), 100));color_range.push_back( std::paircv::Vec3b,double(cv::Vec3b( 250, 180, 140), 200));// define a minimum elevationdouble minElevation -10;// iterate over each pixel in the image, computing the dem pointfor( int y0; yimage.rows; y ){for( int x0; ximage.cols; x ){// convert the pixel coordinate to lat/lon coordinatescv::Point2d coordinate pixel2world( x, y, image.size() );// compute the dem image pixel coordinate from lat/loncv::Point2d dem_coordinate world2dem( coordinate, dem.size() );// extract the elevationdouble dz;if( dem_coordinate.x 0 dem_coordinate.y 0 dem_coordinate.x dem.cols dem_coordinate.y dem.rows ){dz dem.atshort(dem_coordinate);}else{dz minElevation;}// write the pixel value to the fileoutput_dem_flood.atcv::Vec3b(y,x) image.atcv::Vec3b(y,x);// compute the color for the heat map outputcv::Vec3b actualColor get_dem_color(dz);output_dem.atcv::Vec3b(y,x) actualColor;// show effect of a 10 meter increase in ocean levelsif( dz 10 ){add_color( output_dem_flood.atcv::Vec3b(y,x), 90, 0, 0 );}// show effect of a 50 meter increase in ocean levelselse if( dz 50 ){add_color( output_dem_flood.atcv::Vec3b(y,x), 0, 90, 0 );}// show effect of a 100 meter increase in ocean levelselse if( dz 100 ){add_color( output_dem_flood.atcv::Vec3b(y,x), 0, 0, 90 );}}}// print our heat mapcv::imwrite( heat-map.jpg , output_dem );// print the flooding effect imagecv::imwrite( flooded.jpg, output_dem_flood);return 0;
}
解释
在提供的代码片段中有几个关键函数它们负责执行特定的任务如插值、颜色映射和坐标转换。下面是每个关键函数的代码片段和详细解释
1. **线性插值函数 lerp**
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));
} 这个函数执行二维空间中的线性插值。它接受两个点 p1 和 p2以及一个插值比率 t。当 t 从 0 变到 1 时这个函数会生成从 p1 到 p2 的一系列点。
2. **颜色插值模板函数 lerp**
template typename DATATYPE, int N
cv::VecDATATYPE,N lerp(cv::VecDATATYPE,N const minColor,cv::VecDATATYPE,N const maxColor,double const t){cv::VecDATATYPE,N output;for(int i 0; i N; i){output[i] static_castuchar(((1 - t) * minColor[i]) (t * maxColor[i]));}return output;
} 这是一个模板函数用于在两种颜色之间进行插值。它接受两种颜色 minColor 和 maxColor以及插值比率 t。对于颜色中的每个通道对于 cv::Vec3b 是三个通道它计算插值后的颜色值。
3. **计算 DEM 颜色 get_dem_color**
cv::Vec3b get_dem_color(const double elevation){// ... 省略了检查颜色范围和插值的代码 ...return lerp(color_range[idx].first, color_range[idx1].first, t);
} 这个函数根据给定的海拔高度 elevation 从预定义的颜色范围内获取对应的颜色。它首先检查海拔高度是否在颜色范围的边界之外然后找到合适的颜色对进行插值计算。
4. **坐标转换函数 world2dem**
cv::Point2d world2dem(const cv::Point2d coordinate, const cv::Size dem_size){// ... 省略了计算比例和输出点坐标的代码 ...return output;
} 此函数将世界坐标转换为 DEM 图像的像素坐标。它接受一个世界坐标 coordinate 和 DEM 图像的尺寸 dem_size然后计算对应的像素坐标。
5. **坐标转换函数 pixel2world**
cv::Point2d pixel2world(const int x, const int y, const cv::Size size){// ... 省略了计算比例和插值坐标的代码 ...return lerp(leftSide, rightSide, rx);
} 这个函数将像素坐标转换为世界坐标。它接受一个像素坐标 (x, y) 和图像的尺寸 size然后计算出对应的世界坐标。
6. **添加颜色到像素 add_color**
void add_color(cv::Vec3b pix, const uchar b, const uchar g, const uchar r){// ... 省略了颜色添加的代码 ...
} 此函数将特定的颜色值蓝、绿、红通道添加到像素中。它确保添加的颜色值不会超过 255因为颜色值是以 0 到 255 的整数表示的。
7. **主函数 main**
int main(int argc, char* argv[]){// ... 省略了加载图像、DEM 数据和颜色范围设置的代码 ...for(int y 0; y image.rows; y){for(int x 0; x image.cols; x){// ... 省略了像素处理的代码 ...}}// ... 省略了保存图像的代码 ...return 0;
} main 函数是程序的入口点。它首先检查命令行参数然后加载所需的图像和 DEM 数据。接着它通过两个嵌套循环遍历图像的每个像素使用上述函数来计算世界坐标、DEM 坐标、海拔高度和颜色并根据这些信息生成热图和洪水效果图像。
这些函数共同工作实现了一个地理空间数据可视化的程序它可以根据 DEM 数据生成热图并且模拟不同海拔高度下的洪水效果。 如何使用 GDAL 读取栅格数据
此演示使用默认的 OpenCV imread 函数。主要区别在于为了强制 GDAL 加载映像您必须使用适当的标志。 cv::Mat image cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );
加载数字高程模型时每个像素的实际数值是必不可少的不能缩放或截断。例如对于图像数据表示为值为 1 的双精度值的像素与表示为值为 255 的无符号字符的像素具有相同的外观。对于地形数据像素值表示以米为单位的高程。为了确保 OpenCV 保留本机值请在 imread 中使用 GDAL 标志和 ANYDEPTH 标志。 // load the dem modelcv::Mat dem cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );
如果您事先知道要加载的 DEM 模型的类型那么使用断言或其他机制测试 Mattype 或 Mat:d epth 可能是一个安全的选择。NASA 或 DOD 规范文档可以提供各种高程模型的输入类型。主要类型SRTM 和 DTED都是签名短裤。 注意
通常应避免使用经度/纬度地理坐标
地理坐标系是一个球面坐标系这意味着将它们与笛卡尔数学一起使用在技术上是不正确的。此演示使用它们来增加可读性并且足够准确以说明重点。更好的坐标系是通用横轴墨卡托坐标系。
查找拐角坐标
查找图像角坐标的一种简单方法是使用命令行工具 gdalinfo。对于正射校正且包含投影信息的影像可以使用 USGS EarthExplorer。
\f$ gdalinfo N37W123.hgtDriver: SRTMHGT/SRTMHGT File FormatFiles: N37W123.hgtSize is 3601, 3601Coordinate System is:GEOGCS[WGS 84,DATUM[WGS_1984,... more output ...Corner Coordinates:Upper Left (-123.0001389, 38.0001389) (123d 0 0.50W, 38d 0 0.50N)Lower Left (-123.0001389, 36.9998611) (123d 0 0.50W, 36d5959.50N)Upper Right (-121.9998611, 38.0001389) (121d5959.50W, 38d 0 0.50N)Lower Right (-121.9998611, 36.9998611) (121d5959.50W, 36d5959.50N)Center (-122.5000000, 37.5000000) (122d30 0.00W, 37d30 0.00N)... more output ... 结果
以下是程序的输出。使用第一个图像作为输入。对于 DEM 模型请在此处下载位于 USGS 的 SRTM 文件。
http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip 输入图像 热图 热图叠加 参考文献
1、《Reading Geospatial Raster files with GDAL》-----Marvin Smith
