上一节简单介绍了GDAL，这一节将介绍一些GDAL的基本操作，如影像读写、波段提取、波段合成等。代码均用python编写。
1.遥感影像的读写
众所周知，遥感影像是以栅格形式存储的，GDAL中使用dataset表示一个栅格数据（使用抽象类GDALDataset表示），一个dataset包含了对于栅格数据的波段，空间参考以及元数据等信息。一张GeoTIFF遥感影像，一张DEM影像，或者一张土地利用图，在GDAL中都是一个GDALDataset。
- 坐标系统（使用OGC WKT格式表示的空间坐标系统或者投影系统）
- 地理放射变换（使用放射变换表示图上坐标和地理坐标的关系）
- GCPs（大地控制点记录了图上点及其大地坐标的关系，通过多个大地控制点可以重建图上坐标和地理坐标的关系）
- 元数据（键值对的集合，用于记录和影像相关的元数据信息）
- 栅格波段（使用GDALRasterBand类表示，真正用于存储影像栅格值，一个栅格数据可以有多个波段）
- 颜色表（Color Table用于图像显示）
GDAL使用Open（）函数读取影像数据，函数返回Dataset对象。GDAL目前支持约100种格式的栅格数据读取，包括ERDAS Imagine、ENVI、GRASS、GeoTIFF、HDF4、HDF5、TIFF、JPEG、JPEG2000、PNG、GIF、BMP等。
import numpy as np from osgeo import gdal in_ds = gdal.Open(r"D:\data\test.tif") # 打开样本文件 xsize = in_ds.RasterXSize # 获取行列数 ysize = in_ds.RasterYSize bands = in_ds.RasterCount # 获取波段数 geotransform = in_ds.GetGeoTransform() # 获取投影信息 projection = in_ds.GetProjectionRef() block_data = in_ds.ReadAsArray(0,0,xsize,ysize).astype(np.float32)# 获取影像信息 B = block_data[0, :, :] #获取第一波段数据 G = block_data[1,:, :] #获取第二波段数据 R = block_data[2,:, :] #获取第三波段数据
Dataset对象的RasterYSize、RasterXSize和RasterCount属性分别返回栅格数据的行数、列数和波段数。Dataset对象的GetGeoTransform()方法返回栅格数据的坐标转换参数，即行列坐标与空间坐标的转换参数。Dataset对象的GetProjection()方法返回栅格数据的空间参照系统信息（WKT文本）。栅格数据通常有多个波段， Dataset对象的GetRasterBand(index)方法将返回某个波段的数据对象（Band对象），如ds.GetRasterBand(1)返回第一个波段的数据对象（注意：第一波段是1，不是0GetMinimum()和GetMaximum() 方法返回波段数据的最小值和最大值。ComputeStatistics()方法返回波段数据的统计信息，包括最小值、最大值、平均值和标准偏差，approx_ok参数表示是否抽样统计，值为False表示统计所有栅格。ReadAsArray()函数返回的二维数组，astype()可以改变数组的数据类型。
driver = gdal.GetDriverByName('GTiff') out_dataset=driver.Create("RGB.tif",xsize,ysize,3,gdal.GDT_Float32) out_band1 = out_dataset.GetRasterBand(1) out_band1.WriteArray(B) out_band2 = out_dataset.GetRasterBand(2) out_band2.WriteArray(G) out_band3 = out_dataset.GetRasterBand(3) out_band3.WriteArray(R) out_dataset.SetGeoTransform(geotransform) # 写入仿射变换 out_dataset.SetProjection(projection)
将上一步得到的R、G、B三个波段数据存储为新的数据（RGB.tif）。gdal.GDT_Float32代表数据类型，数据类型决定了栅格值的范围，如数据类型为GDT_Byte，则栅格值的范围是0～255；如要存储的数据栅格值的范围是0～65535，则数据类型应该是GDT_UInt16。
gdal常用的数据类型包括：
-
- gdal.GDT_Byte 8bit正整型
- gdal.GDT_UInt16 16bit正整型
- gdal.GDT_Int16 16bit整型
- gdal.GDT_UInt32 32bit 正整型
- gdal.GDT_Int32 32bit整型
- gdal.GDT_Float32 32bit 浮点型
- gdal.GDT_Float64 64bit 浮点型
2.栅格数据行列号和地理坐标相互转换
2.1行列坐标转空间坐标
from osgeo import gdal def Pixel2world(geotransform, line, column): originX = geotransform[0] originY = geotransform[3] pixelWidth = geotransform[1] pixelHeight = geotransform[5] x = column*pixelWidth + originX - pixelWidth/2 y = line*pixelHeight + originY - pixelHeight/2 return(x,y)
line和column为行列坐标（行列号）；pixelWidth和pixelHeigt为x方向比例尺和y方向比例尺；originX和originY为左上角的x和y坐标。
2.2空间坐标转行列坐标：
from osgeo import gdal def world2Pixel(geotransform, x, y): originX = geotransform[0] originY = geotransform[3] pixelWidth = geotransform[1] pixelHeight = geotransform[5] line = int((y-originY)/pixelHeight)+1 column = int((x-originX)/pixelWidth)+1 return (line,column)
3.金字塔生成
图像金子塔是一种多分辨率图像组合的结构。底部是待处理图像的高分辨率表示，顶部是低分辨率的近似。当从下层向上移动时，尺寸和分辨率都会下降。list里面是要生成金字塔的层级，全部为2的幂次，第一层是原始数据。利用gdal生成的ovr数据，实际上就是一个多分辨率的tif文件。文件明在原始的tif后面直接加追.ovr。
# 建立金字塔.ovr ds = gdal.Open(out_path) ds.BuildOverviews('average',[2, 4, 8, 16])
