空间数据处理
July 11, 2018

打开栅格数据的正确方式

打开栅格数据的正确方式

作者：阿振

邮箱：tanzhenyugis@163.com

博客：https://blog.csdn.net/theonegis/article/details/80089375

修改时间：2018-05-16

声明：本文为博主原创文章，转载请注明原文出处

以一个简单例子说明如何打开栅格影像

下面的例子打开一副GeoTIFF影像，输出了影像的一些信息，然后遍历了所有波段，输出波段的一些信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import gdal

# 打开栅格数据集
ds = gdal.Open('example.tif')

# 获得栅格数据的一些重要信息
print(f'投影信息：{ds.GetProjection()}')
print(f'栅格波段数：{ds.RasterCount}')
print(f'栅格列数（宽度）：{ds.RasterXSize}')
print(f'栅格行数（高度）：{ds.RasterYSize}')

# 获取数据集的元数据信息
metadata = ds.GetMetadata_Dict()
for key, value in metadata.items():
    print(f'{key} -> {value}')


for b in range(ds.RasterCount):
    # 注意GDAL中的band计数是从1开始的
    band = ds.GetRasterBand(b + 1)
    # 波段数据的一些信息
    print(f'数据类型：{gdal.GetDataTypeName(band.DataType)}')  # DataType属性返回的是数字
    print(f'NoData值：{band.GetNoDataValue()}')  # 很多影像都是NoData，我们在做数据处理时要特别对待
    print(f'统计值（最大值最小值）：{band.ComputeRasterMinMax()}')  # 有些数据本身就存储了统计信息，有些数据没有需要计算

# 关闭数据集
ds = None

输出如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
投影信息：PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
栅格波段数：3
栅格列数（宽度）：4800
栅格行数（高度）：4800
AREA_OR_POINT -> Area
数据类型：Int16
NoData值：-28672.0
统计值（最大值最小值）：(-435.0, 6134.0)
数据类型：Int16
NoData值：-28672.0
统计值（最大值最小值）：(-468.0, 6265.0)
数据类型：Int16
NoData值：-28672.0
统计值（最大值最小值）：(21.0, 7267.0)

如何将Dataset转为Numpy的ndarray

当我们得到Band对象以后，如果按照GDAL的C/C++接口惯例，我们可以使用WriteRaster()方法进行数据写入（C/C++接口是WriteBlock()），但是在Python中我们有很强大的ndarray对象，所以我们一般是将Band对象中存储的数据转为ndarray进行处理以后，然后再写回去。

下面介绍几种转换的方法：

  1. 在Dataset级别进行转换，转换结果是一个三维数组，第一个维度是波段数
  2. 在Band级别进行转换，转换的结果是一个二维数据
  3. 使用gdal_array模块中的LoadFile()函数直接进行（相当于第一种转换）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import gdal

# 打开栅格数据集
ds = gdal.Open('example.tif')
# 在数据集层面转换
image = ds.ReadAsArray()

print(f'数据的尺寸：{image.shape}')
# 输出结果为：数据的尺寸：(3, 4800, 4800)
# 这说明ReadAsArray方法将每个波段都转换为了一个二维数组

# 获得第一个波段的数据
band1 = image[0]

# 在波段层面的转换
for b in range(ds.RasterCount):
    # 注意GDAL中的band计数是从1开始的
    band = ds.GetRasterBand(b + 1)
    band = band.ReadAsArray()
    print(f'波段大小：{band.shape}')

# 关闭数据集
ds = None

输出结果：

1
2
3
4
数据的尺寸：(3, 4800, 4800)
波段大小：(4800, 4800)
波段大小：(4800, 4800)
波段大小：(4800, 4800)

使用gdal_array模块

1
2
3
4
from osgeo import gdal_array
# gdal_array模块
image = gdal_array.LoadFile('example.tif')
print(f'数据的尺寸：{image.shape}')

在GDAL中使用Python的异常对象

1
2
3
4
5
6
7
8
9
10
11
12
import gdal
import sys

# 允许GDAL跑出Python异常
gdal.UseExceptions()

try:
    ds = gdal.Open('example.tif')
except (FileNotFoundError, RuntimeError) as e:
    print('文件打开失败！')
    print(e)
    sys.exit(1)