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

pip install pillow numpy matplotlib rasterio

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()

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()

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()

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()

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)

# 使用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)

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()}")

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}")