python rasterio包 .tif数据处理 基础操作

最近使用python的rasterio包处理tif数据，记录一下。

官网详细的使用说明：Python Quickstart — rasterio 1.4dev documentation

import rasterio

'''.tif数据基础信息获取'''
tif_path = r"E:/xxx/xx.tif"
img = rasterio.open(tif_path)  # 读取影像
img_bounds = img.bounds  # 影像四至
img_rows, img_cols = img.shape  # 影像行列号
img_bands = img.count  # 影像波段数
img_indexes = img.indexes  # 影像波段
img_crs = img.crs  # 影像坐标系
img_transform = img.transform  # 影像仿射矩阵
upper_left = img.transform * (0, 0)  # 影像左上角像元的左上角点经纬度
lower_right = img.transform * (img.width, img.height)  # 影像右下角像元的右下角点经纬度
img_meta = img.meta  # 影像基础信息，包含driver、数据类型、坐标系、仿射矩阵等

'''.tif数据波段数据获取'''
band1 = img.read(1)  # 读取影像的第一个波段。（波段序号从1开始，不含0波段）
data = img.read([4, 3, 2])  # 读取影像的第4、第3、第2波段。数组的顺序也是432，即data[0]是第4波段，data[1]是第3波段，data[2]是第2波段
data1 = img.read()[2:4, :, :]  # 读取影像的第3、第4波段。不包含第2波段，包含第4波段。
data2 = img.read([4, 3, 2])[:, 150:600, 250:1400]  # 影像切片。读取第150-600行，250-1400列。（含第150行，250列，不含第600行，1400列。数组序号均从0开始）

'''行列号与经纬度间的转换'''
# img.index(lon,lat)获取指定经纬度对应的行列号
# img.xy(row,col)获取指定行列号对应像元的中心经纬度。img.transform * (col, row)获取指定行列号对应像元的左上角经纬度。
# x, y = img.xy(150, 250)，x - img.transform[0] / 2, y - img.transform[4] / 2结果与img.transform * (250, 150)结果一致。
row, col = img.index(img.bounds.left + 1, img.bounds.top - 0.5)  # 获取指定坐标对应的行列号
print(band1[row, col])  # 指定行列号对应的波段数值
x, y = img.xy(150, 250)  # 获取指定行列号对应的像元中心经纬度
center = img.xy(img.height // 2, img.width // 2)  # 影像中心经纬度

'''处理后的数据存储为新的.tif数据'''
# [1]单波段数据存储，数据范围与原始数据一致
new_tif_path = r'D:/test/band1.tif'
new_img = rasterio.open(new_tif_path,
                        'w',
                        driver='GTiff',
                        height=img_rows,
                        width=img_cols,
                        count=1,
                        crs=img.crs,
                        transform=img.transform,
                        dtype=band1.dtype,
                        nodata=0)
new_img.write(band1, 1)  # 写入波段数据
new_img.close()  # 数据存入磁盘

# [2]多波段数据存储，数据范围与原始数据一致
data_path = r'D:/test/data.tif'
with rasterio.open(data_path, 'w', driver='GTiff', height=img_rows, width=img_cols, count=3, crs=img.crs,
                   transform=img.transform, dtype=data.dtype, nodata=0) as dst:
    dst.write(data)  # 或dst.write(data,[1,2,3])
    dst.close()

# [3]切片数据存储，数据范围与原始数据不一致
data2_path = r'D:/test/data2.tif'
out_meta = img.meta.copy()  # 从原始数据中复制影像的基础信息

from rasterio.transform import Affine
transform = Affine.translation(x - img.transform[0] / 2, y - img.transform[4] / 2) * Affine.scale(img.transform[0],
                                                                                                  img.transform[
                                                                                                      4])  # 根据左上角经纬度和分辨率求新的仿射矩阵
out_meta.update({"driver": 'GTiff',
                 "height": data2.shape[1],
                 "width": data2.shape[2],
                 'count': 3,
                 "crs": img.crs,
                 "transform": transform})  # 更新影像信息
with rasterio.open(data2_path, 'w', **out_meta) as dst:
    dst.write(data2)
    dst.close()

原始影像与处理后影像对比：

原始数据                                                                  band1.tif

 

 data.tif                                                                     data2.tif