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

GDAL使用(一)

admin3个月前 (02-03)代码相关119

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

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

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

C++代码:

#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\n", poDataset->GetDriver()->GetDescription(), poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));
	//打印X、Y方向像元数、波段个数
	printf("Size is %d x %d x %d\n", poDataset->GetRasterXSize(), poDataset->GetRasterYSize(), poDataset->GetRasterCount());
	//如果空间参考存在,打印空间参考信息
	if (poDataset->GetProjectionRef() != NULL)
	{
		printf("Projection is `%s'\n", poDataset->GetProjectionRef());
	}
	//获取仿射变换参数并打印
	if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None)
	{
		printf("Origin = (%.6f,%.6f)\n", adfGeoTransform[0], adfGeoTransform[3]);
		printf("Pixel Size = (%.6f,%.6f)\n", 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\n", 
		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\n", 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.\n", pszFormat);
	if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATECOPY, FALSE))
		printf("Driver %s supports CreateCopy() method.\n", pszFormat);


	system("pause");
}


Python代码:

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

@author: cyhu
"""


from osgeo import gdal

# 文件路径
filename = "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))


版权声明:本文由cyhu's essay发布,如需转载请注明出处。

本文链接:https://lovedm.club/?id=115

相关文章

C# 属性,get,set使用

属性反映了状态,是对字段的自然扩展。下面的代码,有一个学生类,类中有年龄属性,通过get,set对年龄进行值的获取与设置,同时对年龄进行约束,使值控制在0-120之间,否则抛出异常信息。namespa...

C# IProgress接口汇报进度

使用IProgress接口来实现进度的汇报,使用CancellationToken类型的参数实现取消操作。IProgress接口中只有一个方法,Report方法,用于报告进度更新。Progress类实...

C# 多线程(2)

C# 多线程(2)

Thread类中的join方法:微软官方解释:Join 一个同步方法,该方法阻止调用线程 (即,调用方法的线程) ,直到 Join 调用方法的线程完成。 使用此...

字节序问题

字节序问题

前阵子读取一个文件,在读取的时候按照文件的描述逐字节读取,但是这个文件的文件头部分数据是分大小端的,这就牵扯到了大小端的转换问题,这里要描述一下字节序。计算机数据存储有两种字节优先顺序:高位字节优先(...

C# WPF程序 显示时间

C# WPF程序 显示时间

例子来源于此视频https://www.bilibili.com/video/av1422127?p=4 此例子是想说明有的类主要使用它的事件namespace Wpf_test { &...

C# 正则表达式(2)

// pattan = @"[^ahou]"; 表示匹配除ahou之外的字符,^在表示反义 string res4 =...