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。

相关文章

使用Docker构建Python Flask程序的详细教程

使用Docker构建Python Flask程序的详细教程

《使用Docker构建PythonFlask程序的详细教程》在当今的软件开发领域,容器化技术正变得越来越流行,而Docker无疑是其中的佼佼者,本文我们就来聊聊如何使用Docker构建一个简单的Py... 目录引言一、准备工作二、创建 Flask 应用程序三、创建 dockerfile四、构建 Docker
阅读更多...
Python使用vllm处理多模态数据的预处理技巧

Python使用vllm处理多模态数据的预处理技巧

《Python使用vllm处理多模态数据的预处理技巧》本文深入探讨了在Python环境下使用vLLM处理多模态数据的预处理技巧,我们将从基础概念出发,详细讲解文本、图像、音频等多模态数据的预处理方法,... 目录1. 背景介绍1.1 目的和范围1.2 预期读者1.3 文档结构概述1.4 术语表1.4.1 核
阅读更多...
Python使用pip工具实现包自动更新的多种方法

Python使用pip工具实现包自动更新的多种方法

《Python使用pip工具实现包自动更新的多种方法》本文深入探讨了使用Python的pip工具实现包自动更新的各种方法和技术,我们将从基础概念开始,逐步介绍手动更新方法、自动化脚本编写、结合CI/C... 目录1. 背景介绍1.1 目的和范围1.2 预期读者1.3 文档结构概述1.4 术语表1.4.1 核
阅读更多...
Conda与Python venv虚拟环境的区别与使用方法详解

Conda与Python venv虚拟环境的区别与使用方法详解

《Conda与Pythonvenv虚拟环境的区别与使用方法详解》随着Python社区的成长,虚拟环境的概念和技术也在不断发展,：本文主要介绍Conda与Pythonvenv虚拟环境的区别与使用... 目录前言一、Conda 与 python venv 的核心区别1. Conda 的特点2. Python v
阅读更多...
Spring Boot中WebSocket常用使用方法详解

Spring Boot中WebSocket常用使用方法详解

《SpringBoot中WebSocket常用使用方法详解》本文从WebSocket的基础概念出发,详细介绍了SpringBoot集成WebSocket的步骤,并重点讲解了常用的使用方法,包括简单消... 目录一、WebSocket基础概念1.1 什么是WebSocket1.2 WebSocket与HTTP
阅读更多...
C#中Guid类使用小结

C#中Guid类使用小结

《C#中Guid类使用小结》本文主要介绍了C#中Guid类用于生成和操作128位的唯一标识符,用于数据库主键及分布式系统,支持通过NewGuid、Parse等方法生成,感兴趣的可以了解一下... 目录前言一、什么是 Guid二、生成 Guid1. 使用 Guid.NewGuid() 方法2. 从字符串创建
阅读更多...
Python使用python-can实现合并BLF文件

Python使用python-can实现合并BLF文件

《Python使用python-can实现合并BLF文件》python-can库是Python生态中专注于CAN总线通信与数据处理的强大工具,本文将使用python-can为BLF文件合并提供高效灵活... 目录一、python-can 库：CAN 数据处理的利器二、BLF 文件合并核心代码解析1. 基础合
阅读更多...
Python使用OpenCV实现获取视频时长的小工具

Python使用OpenCV实现获取视频时长的小工具

《Python使用OpenCV实现获取视频时长的小工具》在处理视频数据时,获取视频的时长是一项常见且基础的需求,本文将详细介绍如何使用Python和OpenCV获取视频时长,并对每一行代码进行深入解析... 目录一、代码实现二、代码解析1. 导入 OpenCV 库2. 定义获取视频时长的函数3. 打开视频文
阅读更多...
Spring IoC 容器的使用详解(最新整理)

Spring IoC 容器的使用详解(最新整理)

《SpringIoC容器的使用详解(最新整理)》文章介绍了Spring框架中的应用分层思想与IoC容器原理,通过分层解耦业务逻辑、数据访问等模块,IoC容器利用@Component注解管理Bean... 目录1. 应用分层2. IoC 的介绍3. IoC 容器的使用3.1. bean 的存储3.2. 方法注
阅读更多...
Python内置函数之classmethod函数使用详解

Python内置函数之classmethod函数使用详解

《Python内置函数之classmethod函数使用详解》：本文主要介绍Python内置函数之classmethod函数使用方式,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地... 目录1. 类方法定义与基本语法2. 类方法 vs 实例方法 vs 静态方法3. 核心特性与用法(1编程客
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2