栅格数据创建与保存
作者:阿振
邮箱:tanzhenyugis@163.com
博客:https://blog.csdn.net/theonegis/article/details/80089375
修改时间:2018-05-24
声明:本文为博主原创文章,转载请注明原文出处
思路与方法
使用Python进行栅格数据处理,很多时候,我们会将GDAL的Dataset
对象转化为NumPy的ndarray
对象,这样我们可以使用很多通用的Python库对数据进行处理,然后再借助GDAL库将数据写回到文件。
不同于普通的二进制文件,空间栅格数据的写需要注意两点:
- 数据的投影信息(确定了平面坐标系)
- 数据的地理坐标信息(确定了图像在给定坐标系下的位置)
在GDAL中,我们首先需要创建Dataset
对象,然后给Dataset
对象填充数据以及元数据信息就OK了。
Driver
或者说GDALDriver
(Python版本的API中对象名称好像都去掉了前缀GDAL,而C/C++版本的API很多对象前面都是有GDAL前缀的,如GDALDataset对象在Python中对应的是Dataset对象)有两个方法:Create()
和CreateCopy()
所以,相应地,我们也有两种思路去创建一个Dataset
对象:
- 如果我们有一个原型数据,比如我们对原始数据进行了处理,处理之后,空间信息,波段等都没有变化,则可以将原始数据作为原型数据,使用
CreateCopy()
方法创建一个和原始数据一样的Dataset
对象,然后在创建好的对象中填充一个ndarray
数据就好了。
- 如果我们没有一个原型数据,那么我们首先需要使用
Create()
方法创建一个空的Dataset
对象,然后手动设置对象的波段,尺寸,空间信息等,然后再在对应的波段填空ndarray
具体的数据。
实现函数
我把上面两种实现思路编码成一个函数,具体实现如下:
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 28 29 30 31 32 33 34 35 36 37 38 39
| 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元数据 """ 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) dataset.SetGeoTransform(transform) dataset.SetProjection(projection) 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的输出数据集。
具体实现如下:
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 28 29 30 31 32 33
| ds = gdal.Open('example.tif')
x_size = ds.RasterXSize y_size = ds.RasterYSize
proj = ds.GetProjection() trans = ds.GetGeoTransform()
image = ds.ReadAsArray()
bnd_red = image[1].astype(float) bnd_nir = image[2].astype(float)
idx_ndvi = (bnd_nir - bnd_red) / (bnd_nir + bnd_red)
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
out2_file = 'DVI.tif'
array2raster(out2_file, idx_ndvi, prototype=gdal.Open(out1_file))
ds = None
|