栅格数据创建与保存

作者：阿振

邮箱：[e-mail]

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

修改时间：2018-05-24

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

思路与方法

使用Python进行栅格数据处理，很多时候，我们会将GDAL的Dataset对象转化为NumPy的ndarray对象，这样我们可以使用很多通用的Python库对数据进行处理，然后再借助GDAL库将数据写回到文件。

不同于普通的二进制文件，空间栅格数据的写需要注意两点：

1. 数据的投影信息（确定了平面坐标系）

2. 数据的地理坐标信息（确定了图像在给定坐标系下的位置）

在GDAL中，我们首先需要创建Dataset对象，然后给Dataset对象填充数据以及元数据信息就OK了。

Driver或者说GDALDriver（Python版本的API中对象名称好像都去掉了前缀GDAL，而C/C++版本的API很多对象前面都是有GDAL前缀的，如GDALDataset对象在Python中对应的是Dataset对象）有两个方法：Create()和CreateCopy()

所以，相应地，我们也有两种思路去创建一个Dataset对象：

1. 如果我们有一个原型数据，比如我们对原始数据进行了处理，处理之后，空间信息，波段等都没有变化，则可以将原始数据作为原型数据，使用CreateCopy()方法创建一个和原始数据一样的Dataset对象，然后在创建好的对象中填充一个ndarray数据就好了。

2. 如果我们没有一个原型数据，那么我们首先需要使用Create()方法创建一个空的Dataset对象，然后手动设置对象的波段，尺寸，空间信息等，然后再在对应的波段填空ndarray具体的数据。

实现函数

我把上面两种实现思路编码成一个函数，具体实现如下：

def array2raster(f_name, np_array, driver='GTiff',
                 prototype=None,
                 xsize=None, ysize=None,
                 transform=None, projection=None,
                 dtype=None, nodata=None):
    """
    将ndarray数组写入到文件中
    :param f_name: 文件路径
    :param np_array: ndarray数组
    :param driver: 文件格式驱动
    :param prototype: 文件原型
    :param xsize: 图像的列数
    :param ysize: 图像的行数
    :param transform: GDAL中的空间转换六参数
    :param projection: 数据的投影信息
    :param dtype: 数据存储的类型
    :param nodata: NoData元数据
    """
    # 创建要写入的数据集（这里假设只有一个波段）
    # 分两种情况：一种给定了数据原型，一种没有给定，需要手动指定Transform和Projection
    driver = gdal.GetDriverByName(driver)
    if prototype:
        dataset = driver.CreateCopy(f_name, prototype)
    else:
        if dtype is None:
            dtype = gdal.GDT_Float32
        if xsize is None:
            xsize = np_array.shape[-1]  # 数组的列数
        if ysize is None:
            ysize = np_array.shape[-2]  # 数组的行数
        dataset = driver.Create(f_name, xsize, ysize, 1, dtype)  # 这里的1指的是一个波段
        dataset.SetGeoTransform(transform)
        dataset.SetProjection(projection)
    # 将array写入文件
    dataset.GetRasterBand(1).WriteArray(np_array)
    if nodata is not None:
        dataset.GetRasterBand(1).SetNoDataValue(nodata)
    dataset.FlushCache()
    return f_name

在使用该函数的时候，要么传进去一个prototype原型数据集，要么传进去transform和projection等信息，这样写入的文件才具有空间参考。

测试案例

下面是一个计算NDVI（Normalized Difference Vegetation Index，归一化植被指数）和DVI（Difference Vegetation Index，差值植被指数）的例子。我们首先计算NDVI，然后通过从原始数据中读取的空间投影和空间变换六元组信息创建输出文件；然后再计算DVI，通过NDVI文件作为原型数据集，以创建DVI的输出数据集。

具体实现如下：

# 打开栅格数据集
ds = gdal.Open('example.tif') # example.tif有三个波段，分别是蓝，红，近红外

# 获取数据集的一些信息
x_size = ds.RasterXSize  # 图像列数
y_size = ds.RasterYSize  # 图像行数

proj = ds.GetProjection()  # 返回的是WKT格式的字符串
trans = ds.GetGeoTransform()  # 返回的是六个参数的tuple

# 在数据集层面ReadAsArray方法将每个波段都转换为了一个二维数组
image = ds.ReadAsArray()

# 获得波段对应的array
bnd_red = image[1].astype(float)  # 红波段
bnd_nir = image[2].astype(float)  # 近红外波段

idx_ndvi = (bnd_nir - bnd_red) / (bnd_nir + bnd_red)  # 计算NDVI指数

out1_file = 'NDVI.tif'
array2raster(out1_file, idx_ndvi,
             xsize=x_size, ysize=y_size,
             transform=trans, projection=proj,
             dtype=gdal.GDT_Float32)

idx_dvi = bnd_nir - bnd_red  # 计算DVI指数

out2_file = 'DVI.tif'
# 这里我们使用out1_file作为原型图像作为参考来保存out2_file
array2raster(out2_file, idx_ndvi, prototype=gdal.Open(out1_file))

# 关闭数据集
ds = None