Python读取和处理分析tif数据的超详细教程_Python

阶段一：基础篇 - 读取和查看tiff数据

1. 安装必要的库

首先安装处理tiff文件常用的库：

pip install pillow numpy matplotlib rasterio

2. 使用pillow读取tiff文件

from pil import image
import matplotlib.pyplot as plt

# 打开tiff文件
img = image.open('example.tif')

# 显示基本信息
print(f"格式: {img.format}")
print(f"大小: {img.size}")
print(f"模式: {img.mode}")

# 显示图像
plt.imshow(img)
plt.show()

3. 转换为numpy数组进行处理

import numpy as np

# 将tiff图像转换为numpy数组
img_array = np.array(img)

# 打印数组信息
print(f"数组形状: {img_array.shape}")
print(f"数据类型: {img_array.dtype}")
print(f"最小值: {img_array.min()}, 最大值: {img_array.max()}")

# 简单处理 - 反转图像
inverted_img = 255 - img_array  # 假设是8位图像
plt.imshow(inverted_img, cmap='gray')
plt.show()

阶段二：中级篇 - 使用rasterio处理地理tiff

1. 读取地理tiff文件

import rasterio

# 打开地理tiff文件
with rasterio.open('geospatial.tif') as src:
    # 读取所有波段数据
    data = src.read()
    
    # 显示元数据
    print(f"波段数: {src.count}")
    print(f"宽度: {src.width}, 高度: {src.height}")
    print(f"坐标参考系统(crs): {src.crs}")
    print(f"地理变换: {src.transform}")
    
    # 显示第一个波段
    plt.imshow(data[0], cmap='viridis')
    plt.colorbar()
    plt.show()

2. 多波段处理和统计

with rasterio.open('multiband.tif') as src:
    # 读取各个波段
    red = src.read(1)   # 红波段
    green = src.read(2) # 绿波段
    blue = src.read(3)  # 蓝波段
    
    # 计算ndvi (假设第4波段是近红外)
    nir = src.read(4)
    ndvi = (nir - red) / (nir + red + 1e-10)  # 避免除以0
    
    # 显示ndvi
    plt.imshow(ndvi, cmap='rdylgn', vmin=-1, vmax=1)
    plt.colorbar(label='ndvi')
    plt.title('ndvi map')
    plt.show()

3. 裁剪和重投影

from rasterio.warp import calculate_default_transform, reproject, resampling

# 裁剪示例
with rasterio.open('large_image.tif') as src:
    # 定义裁剪窗口 (左上x, 左上y, 右下x, 右下y)
    window = rasterio.windows.window(1000, 1000, 500, 500)
    clipped = src.read(window=window)
    
    # 显示裁剪结果
    plt.imshow(clipped[0])
    plt.show()

# 重投影示例
dst_crs = 'epsg:3857'  # web墨卡托投影
with rasterio.open('input.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })
    
    with rasterio.open('reprojected.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=resampling.nearest)

阶段三：高级篇 - 大型tiff处理和优化

1. 分块处理大型tiff

# 使用rasterio的分块读取功能
with rasterio.open('very_large.tif') as src:
    # 获取推荐的块大小
    block_shapes = src.block_shapes
    print(f"块大小: {block_shapes}")
    
    # 迭代处理每个块
    for ji, window in src.block_windows(1):
        # 读取当前块
        block = src.read(window=window)
        
        # 处理块数据 (示例: 计算平均值)
        block_mean = np.mean(block)
        print(f"块 {ji} 的平均值: {block_mean}")
        
        # 这里可以添加你的处理逻辑
        # processed_block = some_processing(block)
        
        # 写入输出文件 (如果需要)
        # dst.write(processed_block, window=window)

2. 使用dask处理超大型tiff

import dask.array as da
import rasterio
from rasterio.windows import window

# 创建延迟加载的dask数组
with rasterio.open('huge.tif') as src:
    # 定义块大小 (例如 1024x1024)
    chunk_size = (1, 1024, 1024)  # (波段, 高, 宽)
    
    # 创建dask数组
    dask_array = da.from_array(src, chunks=chunk_size)
    
    # 现在可以对dask数组执行各种操作 (延迟执行)
    mean_computation = dask_array.mean()
    
    # 实际计算结果 (触发计算)
    print(f"整个图像的平均值: {mean_computation.compute()}")

3. 并行处理多幅tiff

from concurrent.futures import threadpoolexecutor
import glob

def process_tif(file_path):
    """处理单个tiff文件的函数"""
    with rasterio.open(file_path) as src:
        data = src.read(1)
        # 这里添加你的处理逻辑
        result = np.mean(data)
        return file_path, result

# 获取所有tiff文件
tif_files = glob.glob('data/*.tif')

# 使用线程池并行处理
with threadpoolexecutor(max_workers=4) as executor:
    results = list(executor.map(process_tif, tif_files))

# 打印结果
for file_path, mean_val in results:
    print(f"{file_path}: {mean_val}")

实用技巧和小贴士

内存管理:
- 对于大型tiff，使用rasterio的窗口读取或dask进行分块处理
- 处理完成后及时关闭文件或使用with语句

元数据保存:

with rasterio.open('input.tif') as src:
    profile = src.profile  # 保存元数据

# 修改需要的参数
profile.update(dtype=rasterio.float32, count=1)

# 写入新文件时使用保存的元数据
with rasterio.open('output.tif', 'w', **profile) as dst:
    dst.write(processed_data, 1)

tiff格式选择:
- 地理空间数据通常使用geotiff
- 科学数据可能使用bigtiff (用于>4gb的文件)
- 考虑压缩选项 (deflate, lzw等) 以减小文件大小
性能优化:
- 使用rasterio的window参数只读取需要的区域
- 对于多波段操作，考虑使用rasterio的内存映射功能
- 对于重复读取相同文件，考虑转换为更适合分析的格式 (如zarr)