GDAL 操作 TIFF 文件 Python 代码整理：读写、获取坐标系、获取指定位置像元值等教程

GDAL 是一个开源的操作栅格数据和矢量数据的库，本文记录下用 Python 中 GDAL 库操作 TIFF （GeoTIFF）的常见代码，包括读写、获取坐标系、获取指定位置像元值等。

一、TIFF 和 GeoTIFF 介绍

TIFF 简单理解就是一种图像格式，类似于 jpg、png 等。

GeoTIFF 就是在普通 TIFF 文件上增加了地理位置、投影信息、坐标信息等，常用于遥感数据。

二、GDAL 安装

Python 的 GDAL 包我们可以通用 pip 安装，之前老王也分享了安装过程中可能遇到的问题：《pip 安装 GDAL 包出现错误 ERROR: Command errored out with exit status 1 的解决方法》

三、代码分享

以下代码是基于处理 MAIAC AOD 的，所以只有一个波段，因此 Band 数默认就是 1 了，如果是需要多波段，就是 GetRasterBand(band) 的参数需要改一下。

读取 TIFF 文件：

def get_tif_info(tif_path):
    if tif_path.endswith('.tif') or tif_path.endswith('.TIF'):
        dataset = gdal.Open(tif_path)
        pcs = osr.SpatialReference()
        pcs.ImportFromWkt(dataset.GetProjection())
        gcs = pcs.CloneGeogCS()
        extend = dataset.GetGeoTransform()
        # im_width = dataset.RasterXSize #栅格矩阵的列数
        # im_height = dataset.RasterYSize #栅格矩阵的行数
        shape = (dataset.RasterYSize, dataset.RasterXSize)
    else:
        raise "Unsupported file format"

    img = dataset.GetRasterBand(1).ReadAsArray()  # (height, width)
    # img(ndarray), gdal数据集、地理空间坐标系、投影坐标系、栅格影像大小
    return img, dataset, gcs, pcs, extend, shape

经纬度、行列号、投影坐标之间的转换（参考文章里 rowcol_to_xy 写错了）：

def longlat_to_xy(gcs, pcs, lon, lat):
    ct = osr.CoordinateTransformation(gcs, pcs)
    coordinates = ct.TransformPoint(lon, lat)
    return coordinates[0], coordinates[1], coordinates[2]


def xy_to_lonlat(gcs, pcs, x, y):
    ct = osr.CoordinateTransformation(gcs, pcs)
    lon, lat, _ = ct.TransformPoint(x, y)
    return lon, lat


def xy_to_rowcol(extend, x, y):
    a = np.array([[extend[1], extend[2]], [extend[4], extend[5]]])
    b = np.array([x - extend[0], y - extend[3]])

    row_col = np.linalg.solve(a, b)
    row = int(np.floor(row_col[1]))
    col = int(np.floor(row_col[0]))

    return row, col


def rowcol_to_xy(extend, row, col):
    x = extend[0] + col * extend[1] + row * extend[2]
    y = extend[3] + col * extend[4] + row * extend[5]
    return x, y

获取指定位置的值：

def get_value_by_coordinates(tif_pah, coordinates, coordinate_type='rowcol'):
    img, dataset, gcs, pcs, extend, shape = get_tif_info(tif_pah)

    if coordinate_type == 'rowcol':
        value = img[coordinates[0], coordinates[1]]
    elif coordinate_type == 'lonlat':
        x, y, _ = longlat_to_xy(gcs, pcs, coordinates[0], coordinates[1])
        row, col = xy_to_rowcol(extend, x, y)
        value = img[row, col]
    elif coordinate_type == 'xy':
        row, col = xy_to_rowcol(extend, coordinates[0], coordinates[1])
        value = img[row, col]
    else:
        raise 'coordinated_type error'
    return value

保存 TIFF 文件：

def save_tif(array, path, shape, geo_trans, projection, type_name):
    if 'int8' in type_name:
        datatype = gdal.GDT_Byte
    elif 'int16' in type_name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(path, shape[1], shape[0], 1, datatype)

    if dataset is not None:
        dataset.SetGeoTransform(geo_trans)
        dataset.SetProjection(projection)
        dataset.GetRasterBand(1).WriteArray(array)

例子：

1、取值

value = get_value_by_coordinates(tif_path, [115.8597501, 40.0485068], coordinate_type='lonlat')
print(value)

2、保存 TIFF 文件

img_beijing, dataset, gcs, pcs, extend, shape = get_tif_info(aoi_path)

out_path = 'out'

# ratio 是一个 ndarray
save_tif(ratio, out_path + '/ratio.tif', shape, extend, pcs.ExportToWkt(),
         ratio.dtype.name)

参考：

1. https://cloud.tencent.com/developer/news/276097
2. http://zhaoxuhui.top/blog/2019/07/17/GeoTIFFReading&WritingWithGDAL.html