打开栅格数据的正确方式

r囧r小猫 2022-05-24 00:27 282阅读 0赞

# 打开栅格数据的正确方式 #

作者：阿振

邮箱：tanzhenyugis@163.com

博客：[https://blog.csdn.net/theonegis/article/details/80089375][https_blog.csdn.net_theonegis_article_details_80089375]

修改时间：2018-05-16

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

--------------------

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

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

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

输出如下：

投影信息：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()`函数直接进行（相当于第一种转换）

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

输出结果：

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

使用`gdal_array`模块

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

## 在GDAL中使用Python的异常对象 ##

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)

[https_blog.csdn.net_theonegis_article_details_80089375]: https://blog.csdn.net/theonegis/article/details/80089375