当前位置：首页 > 代码相关 > 正文内容

GDAL使用（一）

admin5年前 (2021-02-03)代码相关16952

先开个坑，关于配置什么的有空再写。

这里的代码大多参考官网的教程，可能略有改动。

下面的代码有打开数据、读取栅格相关信息、读取栅格数据、判断栅格驱动是否支持Create()和CreateCopy()等功能。

C++代码：

```cpp
#include <iostream>
#include "gdal_priv.h"
#include "cpl_conv.h" 
#include <cpl_string.h>

int main()
{
	//注册所有驱动
	GDALAllRegister();
	//文件路径
	const char* pszFilename = "E:\\testdata\\landsatdata.tif";
	//打开数据
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(pszFilename, GA_ReadOnly);
	//打开失败返回-1
	if (poDataset == NULL)
	{
		return -1;
	}
	//定义图像仿射变换参数的矩阵
	//六个参数分别为：
	//adfGeoTransform[0] 和 adfGeoTransform[3]分别是图像左上角的x，y坐标
	//adfGeoTransform[1] 和 adfGeoTransform{5]分别是水平和垂直方向上的分辨率，通常两者相等，且垂直方向上的值为负值
	//adfGeoTransform[2] 和 adfGeoTransform[4]表示图像的旋转系数
	double adfGeoTransform[6];
	//打印栅格驱动
	printf("Driver: %s/%s", poDataset->GetDriver()->GetDescription(), poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));
	//打印X、Y方向像元数、波段个数
	printf("Size is %d x %d x %d", poDataset->GetRasterXSize(), poDataset->GetRasterYSize(), poDataset->GetRasterCount());
	//如果空间参考存在，打印空间参考信息
	if (poDataset->GetProjectionRef() != NULL)
	{
		printf("Projection is %s", poDataset->GetProjectionRef());
	}
	//获取仿射变换参数并打印
	if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None)
	{
		printf("Origin = (%.6f,%.6f)", adfGeoTransform[0], adfGeoTransform[3]);
		printf("Pixel Size = (%.6f,%.6f)", adfGeoTransform[1], adfGeoTransform[5]);
	}

int nBlockXSize, nBlockYSize;
	int bGotMin, bGotMax;
	double adfMinMax[2];
	//获取第一个波段信息
	GDALRasterBand* poBand = poDataset->GetRasterBand(1);
	//获取影像的块大小
	poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
	//打印块大小与像元数据类型以及该波段的色彩信息
	printf("Block=%dx%d Type=%s, ColorInterp=%s", 
		nBlockXSize, nBlockYSize, 
		GDALGetDataTypeName(poBand->GetRasterDataType()), 
		GDALGetColorInterpretationName(poBand->GetColorInterpretation()));
	//获取该波段的最大值和最小值
	adfMinMax[0] = poBand->GetMinimum(&bGotMin);
	adfMinMax[1] = poBand->GetMaximum(&bGotMax);
	//printf("%f,%f",adfMinMax[0],adfMinMax[1]);
	if (!(bGotMin && bGotMax))
		//一起获取最大值和最小值
		GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
	//输出最小值和最大值
	printf("Min=%.3fd, Max=%.3f", adfMinMax[0], adfMinMax[1]);

/* 读取栅格数据 */
	//用于存储第一行的数据
	float* pafScanline;
	//获取此波段的X大小
	int nXSize = poBand->GetXSize();
	//安全的分配空间
	pafScanline = (float*)CPLMalloc(sizeof(float) * nXSize);
	//读取数据到分配的空间中（函数的参数比较多）
	poBand->RasterIO(GF_Read, 0, 0, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);

//释放资源
	GDALClose((GDALDatasetH)poDataset);
	CPLFree(pafScanline);

//判断是否支持Create()和CreateCopy()
	const char* pszFormat = "GTiff";
	GDALDriver* poDriver;
	char** papszMetadata;
	poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
	if (poDriver == NULL)
		exit(1);
	papszMetadata = poDriver->GetMetadata();
	if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
		printf("Driver %s supports Create() method.", pszFormat);
	if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATECOPY, FALSE))
		printf("Driver %s supports CreateCopy() method.", pszFormat);

system("pause");
}
```

Python代码：

```python
# -*- coding: utf-8 -*-
"""
Created on Wed Feb  3 15:28:28 2021

@author: cyhu
"""

from osgeo import gdal

# 文件路径
filename = r"E:\testdata\landsatdata.tif"
# 打开数据
dataset = gdal.Open(filename,gdal.GA_ReadOnly)
if not dataset:
    print("open failed!")
# 获取驱动名称
print("Driver:{}/{}".format(dataset.GetDriver().ShortName, dataset.GetDriver().LongName))
# 获取X、Y方向像元数、波段个数
print("Size is {}*{}*{}".format(dataset.RasterXSize, dataset.RasterYSize,dataset.RasterCount))
# 获取空间参考信息
geotransform = dataset.GetGeoTransform()
if geotransform:
    print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
    print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))

# 获取第一个波段
band = dataset.GetRasterBand(1)
# 获取像素数据类型
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))
# 波段的最大值和最小值
min = band.GetMinimum()
max = band.GetMaximum()
if not min or not max:
    (min,max) = band.ComputeRasterMinMax(True)
print("Min={:.3f}, Max={:.3f}".format(min,max))
if band.GetOverviewCount() > 0:
    print("Band has {} overviews".format(band.GetOverviewCount()))
if band.GetRasterColorTable():
    print("Band has a color table with {} entries".format(band.GetRasterColorTable().GetCount()))

# 读取第一行的数据
scanline = band.ReadRaster(xoff = 0, yoff = 0,
                           xsize = band.XSize, ysize = 1,
                           buf_xsize = band.XSize, buf_ysize = 1,
                           buf_type = gdal.GDT_Float32)
# 导入模块转换格式
import struct
tuple_of_floats = struct.unpack('f' * band.XSize, scanline)

# 判断是否支持Create()和CreateCopy()
fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
metadata = driver.GetMetadata()
if metadata.get(gdal.DCAP_CREATE) == "YES":
    print("Driver {} supports Create() method.".format(fileformat))
if metadata.get(gdal.DCAP_CREATECOPY) == "YES":
    print("Driver {} supports CreateCopy() method.".format(fileformat))
```