OpenCV如何使用 GDAL 读取地理空间栅格文件(72)

2024-05-05 23:36
文章标签 读取 使用 opencv 空间 栅格 gdal 72 地理

本文主要是介绍OpenCV如何使用 GDAL 读取地理空间栅格文件(72)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

返回: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 <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 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 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 rangeif( 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 colorreturn 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 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 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 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::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 elevationdouble minElevation = -10;// iterate over each pixel in the image, computing the dem pointfor( int y=0; y<image.rows; y++ ){for( int x=0; x<image.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.at<short>(dem_coordinate);}else{dz = minElevation;}// write the pixel value to the fileoutput_dem_flood.at<cv::Vec3b>(y,x) = image.at<cv::Vec3b>(y,x);// compute the color for the heat map outputcv::Vec3b actualColor = get_dem_color(dz);output_dem.at<cv::Vec3b>(y,x) = actualColor;// show effect of a 10 meter increase in ocean levelsif( dz < 10 ){add_color( output_dem_flood.at<cv::Vec3b>(y,x), 90, 0, 0 );}// show effect of a 50 meter increase in ocean levelselse 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 levelselse if( dz < 100 ){add_color( output_dem_flood.at<cv::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::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] = static_cast<uchar>(((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[idx+1].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 模型的类型，那么使用断言或其他机制测试 Mat：：type（） 或 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.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

输入图像

热图

热图叠加

参考文献：

1、《Reading Geospatial Raster files with GDAL》-----Marvin Smith

这篇关于OpenCV如何使用 GDAL 读取地理空间栅格文件(72)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/962960。

相关文章

sky-take-out项目中Redis的使用示例详解

sky-take-out项目中Redis的使用示例详解

《sky-take-out项目中Redis的使用示例详解》SpringCache是Spring的缓存抽象层,通过注解简化缓存管理,支持Redis等提供者,适用于方法结果缓存、更新和删除操作,但无法实现... 目录Spring Cache主要特性核心注解1.@Cacheable2.@CachePut3.@Ca
阅读更多...
C#下Newtonsoft.Json的具体使用

C#下Newtonsoft.Json的具体使用

《C#下Newtonsoft.Json的具体使用》Newtonsoft.Json是一个非常流行的C#JSON序列化和反序列化库,它可以方便地将C#对象转换为JSON格式,或者将JSON数据解析为C#对... 目录安装 Newtonsoft.json基本用法1. 序列化 C# 对象为 JSON2. 反序列化
阅读更多...
RabbitMQ 延时队列插件安装与使用示例详解(基于 Delayed Message Plugin)

RabbitMQ 延时队列插件安装与使用示例详解(基于 Delayed Message Plugin)

《RabbitMQ延时队列插件安装与使用示例详解(基于DelayedMessagePlugin)》本文详解RabbitMQ通过安装rabbitmq_delayed_message_exchan... 目录 一、什么是 RabbitMQ 延时队列？ 二、安装前准备✅ RabbitMQ 环境要求 三、安装延时队
阅读更多...
Python ORM神器之SQLAlchemy基本使用完全指南

Python ORM神器之SQLAlchemy基本使用完全指南

《PythonORM神器之SQLAlchemy基本使用完全指南》SQLAlchemy是Python主流ORM框架,通过对象化方式简化数据库操作,支持多数据库,提供引擎、会话、模型等核心组件,实现事务... 目录一、什么是SQLAlchemy？二、安装SQLAlchemy三、核心概念1. Engine（引擎）
阅读更多...
Java Stream 并行流简介、使用与注意事项小结

Java Stream 并行流简介、使用与注意事项小结

《JavaStream并行流简介、使用与注意事项小结》Java8并行流基于StreamAPI,利用多核CPU提升计算密集型任务效率,但需注意线程安全、顺序不确定及线程池管理,可通过自定义线程池与C... 目录1. 并行流简介​特点：​2. 并行流的简单使用​示例：并行流的基本使用​3. 配合自定义线程池​示
阅读更多...
GO语言中函数命名返回值的使用

GO语言中函数命名返回值的使用

《GO语言中函数命名返回值的使用》在Go语言中,函数可以为其返回值指定名称,这被称为命名返回值或命名返回参数,这种特性可以使代码更清晰,特别是在返回多个值时,感兴趣的可以了解一下... 目录基本语法函数命名返回特点代码示例命名特点基本语法func functionName(parameters) (nam
阅读更多...
使用shardingsphere实现mysql数据库分片方式

使用shardingsphere实现mysql数据库分片方式

《使用shardingsphere实现mysql数据库分片方式》本文介绍如何使用ShardingSphere-JDBC在SpringBoot中实现MySQL水平分库,涵盖分片策略、路由算法及零侵入配置... 目录一、ShardingSphere 简介1.1 对比1.2 核心概念1.3 Sharding-Sp
阅读更多...
Java 正则表达式的使用实战案例

Java 正则表达式的使用实战案例

《Java正则表达式的使用实战案例》本文详细介绍了Java正则表达式的使用方法,涵盖语法细节、核心类方法、高级特性及实战案例,本文给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要... 目录一、正则表达式语法详解1. 基础字符匹配2. 字符类（[]定义）3. 量词（控制匹配次数）4. 边
阅读更多...
Python Counter 函数使用案例

Python Counter 函数使用案例

《PythonCounter函数使用案例》Counter是collections模块中的一个类,专门用于对可迭代对象中的元素进行计数,接下来通过本文给大家介绍PythonCounter函数使用案例... 目录一、Counter函数概述二、基本使用案例（一）列表元素计数（二）字符串字符计数（三）元组计数三、C
阅读更多...
使用Spring Cache本地缓存示例代码

使用Spring Cache本地缓存示例代码

《使用SpringCache本地缓存示例代码》缓存是提高应用程序性能的重要手段,通过将频繁访问的数据存储在内存中,可以减少数据库访问次数,从而加速数据读取,：本文主要介绍使用SpringCac... 目录一、Spring Cache简介核心特点：二、基础配置1. 添加依赖2. 启用缓存3. 缓存配置方案方案
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2