Python中rasterio库的实现

pip install rasterio

conda install -c conda-forge rasterio

import rasterio

with rasterio.open('sentinel2.tif') as src:
    # 获取数据矩阵（NumPy数组）
    data = src.read()
    
    # 查看元数据
    print(f"数据集信息：\n{src.meta}")
    print(f"坐标系：{src.crs}")
    print(f"影像尺寸：{src.width}x{src.height}")
    print(f"空间分辨率：{src.res}")
    print(f"地理范围：{src.bounds}")

# 创建新栅格文件
profile = src.profile
profile.update(
    dtype=rasterio.float32,
    count=3,
    compress='lzw'
)

with rasterio.open('output.tif', 'w', **profile) as dst:
    dst.write(data.astype(rasterio.float32))

from rasterio.mask import mask
import geojson

# 按像素范围裁剪
window = Window(100, 200, 800, 600)
clip_data = src.read(window=window)

# 按地理范围裁剪，GeoJSON格式
geometry = {
    "type": "Polygon",
    "coordinates": [
        [
            [116.0, 39.5],
            [116.5, 39.5],
            [116.5, 40.0],
            [116.0, 40.0],
            [116.0, 39.5]
        ]
    ]
}

clip_data, clip_transform = rasterio.mask.mask(src, [geometry], crop=True)

from rasterio.enums import Resampling

# 选择重采样方法：最近邻、双线性、立方卷积等
resampled_data = src.read(
    out_shape=(src.height//2, src.width//2),
    resampling=Resampling.bilinear  # 双线性插值法
)

# 地理坐标 ↔ 像素坐标
row, col = src.index(116.3974, 39.9093)  # 天安门坐标转像素位置
lon, lat = src.xy(500, 500)  # 中心点坐标

# 坐标系统转换
from rasterio.warp import transform
dst_crs = 'EPSG:3857'  # Web墨卡托
transformed = transform(src.crs, dst_crs, [lon], [lat])

# 波段合并
with rasterio.open('RGB.tif', 'w', **profile) as dst:
    for i in range(3):
        dst.write(data[i], i+1)

# 计算NDVI
red = src.read(3).astype(float)
nir = src.read(4).astype(float)
ndvi = (nir - red) / (nir + red + 1e-8)

# 可视化NDVI
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 8))
plt.imshow(ndvi, cmap='RdYlGn', vmin=-1, vmax=1)
plt.colorbar(label='NDVI')
plt.title('植被指数分布', fontsize=16)
plt.axis('off')
plt.savefig('ndvi_map.png', dpi=300, bbox_inches='tight')

# 示例：批量处理Landsat8数据提取水体信息

# 1. 读取热红外波段
with rasterio.open('landsat8_b10.tif') as src:
    thermal_band = src.read(1)

# 2. 计算地表温度（简单的示例，实际需要进行辐射校正）
temperature = (thermal_band - 273.15)  # 假设热红外波段已经是开尔文温度

# 3. 应用阈值分割提取水体（基于温度或植被指数）
water_mask = temperature < 20  # 假设温度低于20°C为水体

# 4. 输出二值化结果
with rasterio.open('water_mask.tif', 'w', **src.profile) as dst:
    dst.write(water_mask.astype(rasterio.uint8), 1)

# 5. 生成统计报告
import numpy as np
water_area = np.sum(water_mask) * src.res[0] * src.res[1]  # 计算水体面积
print(f'水体面积: {water_area} 平方米')

from rasterio import Env

with Env(GDAL_NUM_THREADS=4):
    with rasterio.open('large_image.tif') as src:
        data = src.read()