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

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

一、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
赞(2)
关注我们
未经允许不得转载:老王博客 » GDAL 操作 TIFF 文件 Python 代码整理:读写、获取坐标系、获取指定位置像元值等教程

评论 抢沙发

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址