目录
前言
关于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文件数据。