Python遥感数据处理：使用Rasterio读写栅格数据

Rasterio是基于GDAL的Python栅格数据处理库，提供了更Pythonic的API接口。本文将系统介绍Rasterio的常用操作。

一、安装与环境

pip install rasterio matplotlib numpy

二、读取栅格数据

import rasterio
import numpy as np

with rasterio.open("data/sentinel2.tif") as src:
    print(f"CRS: {src.crs}")
    print(f"分辨率: {src.res}")
    print(f"波段数: {src.count}")
    band1 = src.read(1)  # 读取第一波段

三、写入栅格数据

with rasterio.open("output.tif", "w", driver="GTiff",
                    height=band1.shape[0], width=band1.shape[1],
                    count=1, dtype=band1.dtype,
                    crs=src.crs, transform=src.transform) as dst:
    dst.write(band1, 1)

四、NDVI计算

red = src.read(3).astype(float)
nir = src.read(4).astype(float)
ndvi = (nir - red) / (nir + red)

五、栅格裁剪

from rasterio.mask import mask
import json

with open("aoi.geojson") as f:
    geoms = [json.load(f)["features"][0]["geometry"]]

with rasterio.open("input.tif") as src:
    out_image, out_transform = mask(src, geoms, crop=True)