Python GDAL 地学分析（五）读写栅格数据

distance

本节我们介绍使用GDAL读取栅格数据，并将栅格数据转为数组进行处理。

导入第三方包

from osgeo import gdal
import numpy as np

读取栅格图像

函数输入图像所在的路径，输出图像的投影，仿射矩阵以及图像数组。

def read_img(filename):
        dataset = gdal.Open(filename)
        im_width = dataset.RasterXSize #栅格矩阵的列数
        im_height = dataset.RasterYSize #栅格矩阵的行数
        im_geotrans = dataset.GetGeoTransform() #仿射矩阵
        im_proj = dataset.GetProjection() #地图投影
        im_data = dataset.ReadAsArray(0,0,im_width,im_height)
        del dataset
        return im_proj,im_geotrans,im_data

写入栅格图像

输入要保存的路径，投影，仿射矩阵，图像数组，输出为指定格式的栅格影像。

def write_img(filename,im_proj,im_geotrans,im_data):
        if "int8" in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        if "int16" in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape
        driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影