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 gdal

os.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])：

1. filename：创建的数据集的路径
2. ysize：新数据集的列数
3. bands：新数据集里的波段数，默认值为1
4. data_type：将要储存在新数据集中的数据类型，默认为GDT_Byte
5. 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])

1. xoff：列读取的起点，默认值为0
2. yoff：行读取的起点，默认值为0
3. win_xsize：要读取的列数，默认读取所有列
4. win_ysize：要读取的行数，默认读取所有行
5. buf_xsize：输出数组里的列数，默认为win_xsize的值，如果值不同于win_xsize，数据将会重采样
6. buf_ysize：输出数组里的行数，默认为win_xsize的值，如果值不同于win_ysize，数据将会重采样
7. 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 gdal


data_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 np

data = 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 gdal

os.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_xsize
    else:
        cols = xsize - x

    # 在y方向循环blocks
    for y in range(0, ysize, block_ysize):

        # 获取要读取的行数
        if y + block_ysize < ysize:
            rows = block_ysize
        else:
            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 gdal

data_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）处的像素值为：

62

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

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

62

注意： 使用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, 5267134

os.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

  结果展示：

  图像属性信息：