【实战】基于GDAL库读取指定经纬度下的地表覆盖数据(数据源：清华大学FROM_GLC10(2017))

2023-11-15

前言

关于GDAL库的网上已经有很多版本了，自行下载合适的版本即可。
提醒：GDAL依赖proj和sqlite库。
相关链接：
SQLite3源码下载与编译（开发环境：Win10+VS2022）
PROJ 9.1.1源码下载编译（Win10+VS2022）

数据源

清华大学宫鹏教授学科组10m土地覆盖数据

FROM_GLC10(2017)数据由清华大学宫鹏教授学科组，基于10m分辨率的Sentinel-2影像制作，首个10米分辨率的全球土地覆盖产品。新产品依据全球约9.3万个样本点位上30余万个不同季节的训练样本所训练的分类器，通过对样本数量和误差的的深入分析和模拟，提出“地表覆盖有限样本稳定分类”理论。FROM-GLC10 产品的整体准确率为 72.76%。
在这里插入图片描述

数据集类型

数据集类型为常规的十种，具体见下列枚举。

enum class LandCover
{
	Error = 0,				// 错误
	Cropland = 10,			// 耕地
	Forest = 20,			// 林地
	Grassland = 30,			// 草地
	Shrubland = 40,			// 灌木
	Wetland = 50,			// 湿地
	Water = 60,				// 水域
	Tundra = 70,			// 冻土
	ImperviousSurface = 80, // 不透水面 
	Bareland = 90,			// 裸地
	SnowLce = 100,			// 雪/冰
};

下载途径

由于所有图片需要一个个下载很繁琐，推荐将下载网页保存至本地，然后通过正则提取所有的下载链接，然后通过迅雷批量下载，一次可以下载1000个。总共大约120g。官方下载链接

GDAL库读取FROM_GLC10数据集

关于GDAL读取tiff格式图片数据我前面已经介绍过了。请移步：使用GDAL库读取Tiff文件

下载一个GIS平台

这里我推荐图新地球，简单绿色，解压即可使用。。图新地球 4（LocaSpaceViewer）官网

在这里插入图片描述

加载影像数据只需要将图片拖至图层即可。双击即可实现跳转。

数据集命名规则

通过查看图片的范围以及图片命名可以分析出。
在这里插入图片描述

图片命名格式：fromglc10v01_维度_经度.tif

而且每张图片都是2°的正方形（已经验证成功）。
所以我们可以根据所提供的经纬度计算出该点所在的文件位置。具体公式如下：

const int y = static_cast<int>((lat + 180) / 2) * 2 - 180;
const int x = static_cast<int>((lon + 180) / 2) * 2 - 180;
const std::string fileName = m_path + "/fromglc10v01_" + std::to_string(y) + "_" + std::to_string(x) + ".tif";

GetGeoTransfrom方法介绍

我们可以通过GetGeoTransfrom读取该图片的一些重要坐标信息。

实例代码

#include "gdal_priv.h"
GDALDataset *ds = (GDALDataset*)GDALOpen(filename, GA_ReadOnly);
double geoTransform[6];
ds->GetGeoTransform(geoTransform);

geoTransform参数说明

geoTransform共包括6个参数，具体见下方说明

下标	说明
GeoTransform[0]	左上角横坐标（投影坐标或者是实际的经纬度这取决于数据集）
GeoTransform[1]	像元宽度(影像在水平空间的分辨率)
GeoTransform[2]	行旋转
GeoTransform[3]	左上角纵坐标（投影坐标或者是实际的经纬度这取决于数据集）
GeoTransform[4]	列旋转
GeoTransform[5]	像元高度(影像在垂直空间的分辨率)

说明：如果影像是指北的，GeoTransform[2]和GeoTransform[4]这两个参数的值为0，GeoTransform[5]为负；
如果图像不含地理坐标信息，默认返回值是：(0,1,0,0,0,1)

计算出影像图片的经纬度范围

根据geoTransform方法我们可以得出左上角横坐标以及左上角纵坐标。
以fromglc10v01_0_100.tif图片为例。

输出的信息为

99.9999
8.98315e-05
0
2.00001
0
-8.98315e-05

可以看到刚好为图片的左上角经纬度，由此可得出该数据集存储的信息为实际的经纬度。那么读取就更加简单了。

根据图片的分辨率为22265*22265 ，我们此处证明前面说的图片范围为2°：
图片实际宽度 = 分辨率宽度 * 像元宽度 = 22265 * 8.98315e-05 = 2.01°

在这里插入图片描述

那么我们可以算出整张图片的四个顶点：

	std::array<std::pair<int,int>, 4> pos; // 左上 右上 右下 左下
	pos.at(0) = std::pair<int, int>(static_cast<int>(round(adfGeoTransform[0])), 
		static_cast<int>(round(adfGeoTransform[3])));
	pos.at(1) = std::pair<int, int>(static_cast<int>(round(adfGeoTransform[0])),
		static_cast<int>(round(adfGeoTransform[3] + (m_poDataset->GetRasterXSize() - 1) * adfGeoTransform[5])));
	pos.at(2) = std::pair<int, int>(static_cast<int>(round(adfGeoTransform[0] + (m_poDataset->GetRasterYSize() - 1) * adfGeoTransform[1])),
		static_cast<int>(round(adfGeoTransform[3])));
	pos.at(3) = std::pair<int, int>(static_cast<int>(round(adfGeoTransform[0] + (m_poDataset->GetRasterYSize() - 1) * adfGeoTransform[1] + (m_poDataset->GetRasterXSize() - 1) * adfGeoTransform[2])),
		static_cast<int>(round(adfGeoTransform[3] + (m_poDataset->GetRasterYSize() - 1) * adfGeoTransform[4] + (m_poDataset->GetRasterXSize() - 1) * adfGeoTransform[5])));

根据经纬度定位图中像素位置

此时我们已经知道图中分辨、实际经纬度大小、实际图片经纬度范围，那么公式如下：

图中像素坐标 = 图中分辨率/实际经纬度大小 * 实际图片经纬度范围

逼近算法

if (lon >= pos.at(0).first && lon <= pos.at(2).first && lat <= pos.at(0).second && lat >= pos.at(3).second)
{
		// 计算出来的所在行列
		double x{ .0 };
		double y{ .0 };

		// 计算出所在的经纬度
		double px{ lon };
		double py{ lat };

		// 逼近算法
		do
		{
			y = (pos[2].first - px) / adfGeoTransform[1];  
			x = (pos[1].second - py) / adfGeoTransform[5]; 

			px = adfGeoTransform[0] + x * adfGeoTransform[1] + y * adfGeoTransform[2];
			py = adfGeoTransform[3] + x * adfGeoTransform[4] + y * adfGeoTransform[5];

		} while (fabs(px - lon) > 0.0001 || fabs(py - lat) > 0.0001);

		
}

读取指定像素坐标的像元值

这里我们需要了解GDAL强大的一个方法RasterIO；

函数原型如下：

    CPLErr      RasterIO( GDALRWFlag, int, int, int, int,
                          void *, int, int, GDALDataType,
                          GSpacing, GSpacing, GDALRasterIOExtraArg* psExtraArg
#ifndef DOXYGEN_SKIP
                          OPTIONAL_OUTSIDE_GDAL(nullptr)
#endif
                          ) CPL_WARN_UNUSED_RESULT;

参数介绍

参数	介绍
参数一	读写标记来指定你是读取图像还是写入图像
参数二参数三	读取的位置，起始行列号(像素坐标)
参数四参数五	读写图像的宽度和高度例如本文只需要读取一个像元值就是 1*1
参数六	存储图像的元素值的地方读取图像时表示读取出的像元值存放位置，写入图像时表示该数据将会写入指定位置上去
参数七参数八	表示参数六数组的大小，x*y 还可以放缩图片对原有区域进行重采样。如果这俩个参数大于缓冲区（参数六）那么将会返回错误值
参数九	将数据读取的数据类型为什么，这个取决于参数六的数据类型。例如参数六为int[],那么此处填写GDT_CInt32即可
参数十参数十一	控制参数六中像元的存储顺序，参数十表示的是当前像素值和下一个像素值之间的间隔，参数十一表示当前行和下一行的间隔。单位都是按照字节为单位计算。一般默认给0，0

std::unique_ptr<int> pixedValue(new int[2]);
		CPLErr pd = poBand->RasterIO(GF_Read, x, y, 1, 1,
			pixedValue.get(), 1, 1, GDT_CInt32, 0, 0);

读取指定经纬度的地表覆盖类型

通过查看GIS平台可以看到该图片在左上角(101.9,1.8)位置应该为水域。
在这里插入图片描述

测试结果
60的枚举为水域，符合GIS显示结果。
在这里插入图片描述

最后

本文完，后续介绍如何解析带有投影坐标的tiff文件数据。

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)