Python地理数据处理十二：栅格数据读写

1.GDAL

GDAL（Geospatial Data Abstraction Library ）是非常受欢迎的、强大的栅格文件读写库。GDAL库是开源的，但是有宽松的授权，所以许多商业软件包都使用它。现有的大部分GIS或者遥感平台，不论是商业软件ArcGIS，ENVI还是开源软件GRASS，QGIS，都使用了GDAL作为底层构建库。
GDAL库以其能读写很多不同格式而被人熟知，但是也包含一些数据处理功能，如邻近分析。NumPy模块是专门处理大矩阵数据，并且可以使用GDAL直接将数据读进NumPy数组。在操作数据之后，可以使用NumPy或其他模块来处理这些矩阵，可以将数组写回磁盘作为一个栅格数据集。
可以查看GDAL光栅驱动：

每个驱动程序读写特定的数据格式。GDAL的数据结构大体上和栅格数据集匹配，每个数据集包含一个或多个波段，每个波段又包含像素数据和可能的概视图。地理参考信息包含在数据集中，因为所有的波段都使用相同的地理参考信息。

美国陆地卫星（LANDSAT）系列卫星由美国航空航天局（NASA）和美国地质调查局（USGS）共同管理。自1972年起，LANDSAT系列卫星陆续发射8颗（第6颗发射失败）。Landsat1-5目前均已退役，landsat7则自2003年6月以来，该传感器的扫描线校正器（SLC）出现故障导致的数据间隙数据，但目前仍在运行。Landsat8 于2013年2月11日发射，经过100天测试运行后开始获取影像，目前仍在运。
美国地质调查局将Landsat影像处理成 GeoTIFFs 影像集。除了波段6（热红外）和波段8（全色），每一个波段都是30米的分辨率，因为来自于同一个Landsat场景，具有相同的尺寸。波段可以直接叠置在另一个上面，而不需要修改。

1.将独立的栅格波段合成一张图像：

# Script to stack three 1-band rasters into a a single 3-band image.
import os
from osgeo import gdalos.chdir(r'E:\gis with python\Landsat\Washington')
band1_fn = 'p047r027_7t20000730_z10_nn10.tif'
band2_fn = 'p047r027_7t20000730_z10_nn20.tif'
band3_fn = 'p047r027_7t20000730_z10_nn30.tif'# 打开GeoTIFF波段1，索引为1
# 波段索引从1开始，而不是0
in_ds = gdal.Open(band1_fn)
in_band = in_ds.GetRasterBand(1)# 利用和波段1相同的属性创建三波段GeoTIFF
gtiff_driver = gdal.GetDriverByName('GTiff')  # 使用驱动对象创建新的数据集
out_ds = gtiff_driver.Create('nat_color.tif',in_band.XSize, in_band.YSize, 3, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())# 从输入波段复制像素到输出波段3
in_data = in_band.ReadAsArray()
out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)# 从数据集，而不是波段复制像素
in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())# 为每个波段计算统计值
out_ds.GetRasterBand(1).WriteArray(gdal.Open(band3_fn).ReadAsArray())out_ds.FlushCache()
for i in range(1, 4):out_ds.GetRasterBand(i).ComputeStatistics(False)# 创建概视图或金字塔图层
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])del out_ds

RGB图像：

driver.Create(filename, xsize, ysize, [bands], [data_type], [options])：

filename：创建的数据集的路径
ysize：新数据集的列数
bands：新数据集里的波段数，默认值为1
data_type：将要储存在新数据集中的数据类型，默认为GDT_Byte
option：创建选项字符串列表。值是基于所创建的数据集类型。

GDAL数据类型：

常量	数据类型
GDT_Unknown	Unknown
GDT_Byte	Unsigned 8-bit integer (byte)
GDT_UInt16	Unsigned 16-bit integer
GDT_Int16	Signed 16-bit integer
GDT_UInt32	Unsigned 32-bit integer
GDT_Int32	Signed 32-bit integer
GDT_Float32	32-bit floating point
GDT_Float64	64-bit floating point
GDT_CInt16	16-bit complex integer
GDT_CInt32	32-bit complex integer
GDT_CFloat32	32-bit complex floating point
GDT_CFloat64	64-bit complex floating point
GDT_TypeCount	Number of available data types

从输入波段中获取此信息，但是这些图像使用了GDT_Byte默认类型。

1. 查看空的三波段数据集的SRS：

out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())

从输入数据集中得到SRS并复制到新数据集，再对 geotransform 做同样的操作。geotransform 提供原始坐标和像素大小，并伴随着旋转值，可使图像顶部朝北。如果需要将数据集放置到正确的空间位置，数据集和像素大小是非常重要的参数。

2. 创建数据集后，就要添加像素值，将读取的波段1的像素值读进 NumPy 数组。不给 Readasarray（） 函数提供任何参数，所有的像素值就会返回在一个和栅格具有相同尺寸的二维数组里。in_data变量保存了一个二维数组像素值：

in_data = in_band.ReadAsArray()

3. 因为波段1是蓝色波段，所以必须放到输出图层的第3个波段位置上，从而得到RGB波段顺序。
然后，从out_data获取第3个额伯顿啊，并用WriteArray（）函数将in_data数组里的像素值复制到新数据集的第3波段。

out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)

4. 将红色和绿色波段添加到数据集中，并打开第2波段的GeoTIFF。直接从数据集中读取像素数据。

in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())

当数据集上调用ReadAsArray函数时，如果正在读取的数据集有多个波段，将会得到一个三维数组。因为Landsat文件只有一个波段，所以返回的是二维数组。

5. 对红色像素值做同样的事情，但是压缩为更少的代码：

out_ds.GetRasterBand(1).WriteArray(gdal.Open(band3_fn).ReadAsArray())

6. 计算数据集中每个波段的统计数据，能够使软件更容易显示。统计数据包括波段的平均值、最小值、最大值和标准查。使得GIS软件可以使用这些信息来拉伸屏幕数据，使得看起来更好。
统计数据之前，必须确保数据已经写入磁盘，计算每个波段的统计数据。传递False给这个函数，告诉它你需要的实际的统计数据，而不是估计值，因为他可能从概视图中或采样像素的子集中得到。若估计值可以接受，就传递True给函数；使得计算速度更快，因为不是每个像素都需要检查：

out_ds.FlushCache()
for i in range(1, 4):out_ds.GetRasterBand(i).ComputeStatistics(False)

7. 建立数据集的概视图

out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])

可以在rasterio上查看，此模块依赖于GDAL，实质上是使用GDAL来读写数据，但rasterio模块尝试使栅格数据的处理工作更加轻松。
imageio模块不依赖于GDAL，是纯用python编写的模块，不专注于地理空间数据，可以读取和写入许多不同的栅格格式，包括视频格式。

2.读取部分数据集

访问子集，而不是访问整个图像。ReadAsAsrray（）函数有几个可选参数，会根据你是否在用数据集或波段而变化。
band.ReadAsArray([xoff], [yoff], [win_xsize], [win_ysize], [buf_xsize], [buf_ysize], [buf_obj])

xoff：列读取的起点，默认值为0
yoff：行读取的起点，默认值为0
win_xsize：要读取的列数，默认读取所有列
win_ysize：要读取的行数，默认读取所有行
buf_xsize：输出数组里的列数，默认为win_xsize的值，如果值不同于win_xsize，数据将会重采样
buf_ysize：输出数组里的行数，默认为win_xsize的值，如果值不同于win_ysize，数据将会重采样
buf_obj：用于存放数组，而不是创建数组的NumPy数组。如果需要，数据将会重采样来适应数组。值将会转换为这种数组的类型。

xoff和yoff参数分别用于指定开始读取的列和行的偏移值。默认从第一行和第一列开始。
win_xsize和win_ysize参数用来显示读取了多少行和列，默认读取所有行和列
如果提供的buf_xsize和buf_ysize和该数组的大小不匹配，就会得到一个错误，但是没必要给出buf_xsize和buf_ysize的值，因为这用数组本身确定。

1. 读取从第6000行第1400列开始的3行和6列数据：

import os
import numpy as np
from osgeo import gdaldata_dir = r'E:\gis with python'# 打开Landsat波段
os.chdir(os.path.join(data_dir, 'Landsat', 'Washington'))
ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
band = ds.GetRasterBand(1)# 读取指定范围数据
data = band.ReadAsArray(1400, 6000, 6, 3)
print(data)

2. 数据类型转换：

# 使用numpy将数据转换为浮点数
data = band.ReadAsArray(1400, 6000, 6, 3).astype(float)
print(data)

3. 读取数据时，用GDAL来转换：
首先，创建一个浮点数组，然后把它作为一个buf_obj参数传递给readasarray，同时确保创建的数组于正在读取的数据大小相同。

# 通过将它们读入浮点数组，来转换为浮点数
import numpy as npdata = np.empty((3, 6), dtype=float)  # dtype参数用来指定数组的数据类型
band.ReadAsArray(1400, 6000, 6, 3, buf_obj=data)
print(data)

4. 将数据数组写入其他数据集的指定位置，给函数传递偏移值。它将从提供的偏移值开始，写入传递给 WriteArray（） 函数的数组里的所有数据：

band2.WriteArray(data, 1400, 6000)

5.处理一个超过内存大小的大数据集： 一次只处理一块。
注意： 栅格数据在磁盘上按块储存，所以在这些块中处理图像效率非常高。如图所示：

将DEM的高程从米制转换成英尺（一次只处理一块）：

import os
import numpy as np
from osgeo import gdalos.chdir(r'E:\gis with python\Washington\dem')# 打开输入光栅并获得其尺寸
in_ds = gdal.Open('gt30w140n90.tif')
in_band = in_ds.GetRasterBand(1)
xsize = in_band.XSize
ysize = in_band.YSize# 获取块大小和NoData值
block_xsize, block_ysize = in_band.GetBlockSize()
nodata = in_band.GetNoDataValue()# 创建具有相同维度和数据类型的输出文件
out_ds = in_ds.GetDriver().Create('dem_feet.tif', xsize, ysize, 1, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)# 在x方向循环这些块
for x in range(0, xsize, block_xsize):# 获取要读取的列数if x + block_xsize < xsize:cols = block_xsizeelse:cols = xsize - x# 在y方向循环blocksfor y in range(0, ysize, block_ysize):# 获取要读取的行数if y + block_ysize < ysize:rows = block_ysizeelse:rows = ysize - y# 读取一个块的数据值，将其转换为英尺# 然后将结果写入输出中相同的块位置data = in_band.ReadAsArray(x, y, cols, rows)data = np.where(data == nodata, nodata, data * 3.28084)out_band.WriteArray(data, x, y)# 刷新缓存
# 设置NoData值后进行统计
out_band.FlushCache()
out_band.SetNoDataValue(nodata)
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
del out_ds

这种方法比一次性读取整个波段来说更加复杂，但是节约内存。

3.现实世界坐标

前提：真实世界坐标和栅格数据坐标使用相同的SRS

需要的数据：坐标原点、像素大小（单个像元所代表的实际地物大小）和旋转值，都储存在从数据集间复制过来的 geotransform（地理变换）里。地理变换是一个元组，包含六个值：

索引	描述
0	原点的x坐标
1	像素宽度
2	x像素旋转（图像朝北为0°）（默认为0）
3	原点的y坐标
4	y像素旋转（图像0°为朝北）（默认为0）
5	像素高度（负值）

GDAL提供的ApplyGeoTransform函数来进行仿射变换，需要geotransform、x值和y值3个参数。当使用一个数据集的geotransform变换时，这个函数能将图像坐标（偏移）转换为现实世界坐标。
GDAL提供了一个ApplyGeoTransform函数进行仿射变换，parameter：geotransform、x值、y值。

另一种思路（数据集的逆变换）：

1. GDAL 1.X版本： InvGeoTransform函数

# 现在得到逆变换
# 原始文件可用于将偏移量转换为真实世界的坐标
# 而逆文件可用于将真实世界的坐标转换为偏移量import os
import numpy as np
from osgeo import gdaldata_dir = r'D:\geodata'#从陆地卫星的一个波段中获取geotransform
os.chdir(os.path.join(data_dir, 'Landsat', 'Washington'))
ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
band = ds.GetRasterBand(1)
gt = ds.GetGeoTransform()
print(gt)# GDAL 1.x: 得到一个成功标志和地理变换
success, inv_gt = gdal.InvGeoTransform(gt)
print(success, inv_gt)

成功返回1，不成功返回0。 0：表示仿射变换不可逆。

2. GDAL 2.X版本（我的版本为3.2.1）： InvGeoTransform函数，使用此函数将现实世界坐标转换为图像坐标。

inv_gt = gdal.InvGeoTransform(gt)
print(inv_gt)

(-12060.5, 0.03508771929824561, 0.0, 188406.5, 0.0, -0.03508771929824561)

对照上面的表格进行解释。

3. 坐标转换实例 1： 获取（465200, 5296000）处的像素值

# 使用inverset geotransform从真实世界的UTM坐标获得像素偏移量
# 输出为浮点型
offsets = gdal.ApplyGeoTransform(inv_gt, 465200, 5296000)
print(offsets)# 将偏移量转换为整型
xoff, yoff = map(int, offsets)
print(xoff, yoff)

[4262.307017543859, 2581.9385964912362]
4262 2581

# 读取像素值
# ReadAsArray函数需要使用整型
# ReadAsArray函数返回一个二维数组
value = band.ReadAsArray(xoff, yoff, 1, 1)[0,0] # [0,0]表示数组里的第一行，第一列的值
print(value)

最终，（465200, 5296000）处的像素值为：

存在的问题： 单个像素值的提取效率低，不适用于整个波段的像素值。
改进的方法：

# 先读取所有像素值
# 再取出所需要的像素值
data = band.ReadAsArray()
x, y = map(int, gdal.ApplyGeoTransform(inv_gt, 465200, 5296000))
value = data[yoff, xoff] # 这里需要注意，这里使用的是Numpy读取的
print(value)

注意： 使用Numpy数组的[行，列]偏移与GDAL使用的方式相反。

4. 实例 2： 提取图像空间子集，并保存（感兴趣区域已知）。即将图像的现实世界坐标转换为像素偏移值。
数据说明： 一幅真彩色Landsat图像（Landsat_color.tif）

import os
from osgeo import gdal# 提取的边框的坐标
interested_ulx, interested_uly = 517296, 5296820
interested_lrx, interested_lry = 540315, 5267134os.chdir(r'D:\geodata\Landsat\Washington')
in_ds = gdal.Open('Landsat_color.tif')# 将现实世界的坐标转换为像素偏移量
in_gt = in_ds.GetGeoTransform()
inv_gt = gdal.InvGeoTransform(in_gt)
if gdal.VersionInfo()[0] == '1':if inv_gt[0] == 1:inv_gt = inv_gt[1]else:raise RuntimeError('Inverse geotransform failed')
elif inv_gt is None:raise RuntimeError('Inverse geotransform failed')# 获取与边框角坐标对应的偏移量
# 计算左上角和右下角偏移
offsets_ul = gdal.ApplyGeoTransform(inv_gt, interested_ulx, interested_uly)
offsets_lr = gdal.ApplyGeoTransform(inv_gt, interested_lrx, interested_lry)
off_ulx, off_uly = map(int, offsets_ul)
off_lrx, off_lry = map(int, offsets_lr)# 计算需要提取的行列数
rows = off_lry - off_uly
columns = off_lrx - off_ulx# 创建具有正确行列数的输出光栅
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('interested.tif', columns, rows, 3)
out_ds.SetProjection(in_ds.GetProjection())# 将新的原点坐标放入地理变换中
subset_ulx, subset_uly = gdal.ApplyGeoTransform(in_gt, off_ulx, off_uly)
out_gt = list(in_gt)
out_gt[0] = subset_ulx
out_gt[3] = subset_uly
out_ds.SetGeoTransform(out_gt)# 遍历三个波段
for i in range(1, 4):in_band = in_ds.GetRasterBand(i)out_band = out_ds.GetRasterBand(i)# 利用计算值读入数据data = in_band.ReadAsArray(off_ulx, off_uly, columns, rows)# 从原点处写入数据out_band.WriteArray(data)del out_ds

结果展示：

图像属性信息：