使用Python高效解析东北10米精度作物时序栅格数据
从高分辨率遥感数据到可操作农业洞察：实战指南
在处理覆盖中国东北地区（2017–2024年）的10米空间分辨率作物分类栅格数据时，仅掌握基础读取方法远远不够。真正的挑战在于如何高效、稳健地完成跨年度数据的统一处理与分析。本文聚焦于利用现代Python工具链，构建一套可复用、可扩展的数据处理流水线。
环境配置：构建稳定地理分析环境
为避免依赖冲突，推荐使用Conda创建专用环境：
conda create -n agri_raster python=3.9
conda activate agri_raster
conda install -c conda-forge rasterio geopandas gdal numpy pandas matplotlib jupyterlab
核心库说明：
rasterio：专用于读写栅格数据，支持多种格式及投影转换。geopandas：用于加载行政区划边界等矢量参考数据。numpy：执行高效的数组运算。matplotlib：生成高质量可视化图表。
数据初探：建立元信息认知
打开任意一年的TIF文件，获取关键元数据以评估数据质量：
import rasterio
import numpy as np
with rasterio.open("data/NE_Crop_2023.tif") as src:
print("图像尺寸:", src.shape)
print("地理范围 (左上, 右下):", src.bounds)
print("坐标系统:", src.crs)
print("像素大小:", src.res) # (x_res, y_res)
print("波段类型:", src.dtypes[0])
print("有效值范围:", np.unique(src.read(1)))
输出示例：
- 尺寸:
(120000, 96000) - 坐标系:
EPSG:4326 - 分辨率:
(10.0, 10.0)米 - 有效标签:
[0, 1, 2, 3]→ 对应空地、水稻、玉米、大豆
⚠️ 提示：对于超大文件，避免一次性读取整图。应优先通过
src.profile查看头信息，或使用read(window=(row_start, col_start, height, width))实现分块读取。
统一坐标系统一：确保时空对齐
若不同年份的栅格坐标系不一致，需进行重投影。以下为标准化流程：
from rasterio.warp import reproject, Resampling
from rasterio.transform import Affine
def align_crs(input_path, output_path, target_crs="EPSG:4326"):
with rasterio.open(input_path) as src:
profile = src.profile.copy()
profile.update({
'crs': target_crs,
'transform': Affine.identity(),
'width': src.width,
'height': src.height
})
with rasterio.open(output_path, 'w', **profile) as dst:
for i in range(1, src.count + 1):
data = src.read(i)
reproject(
source=data,
destination=dst.write(data, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=dst.transform,
dst_crs=target_crs,
resampling=Resampling.nearest
)
调用此函数将所有年份数据统一至同一坐标系，为后续时间序列分析打下基础。
时序分析：提取作物面积动态变化
基于统一后的数据，统计每年各类作物的总面积（单位：平方米），并导出为时间序列表格：
import pandas as pd
years = range(2017, 2025)
results = []
for year in years:
path = f"data/NE_Crop_{year}.tif"
with rasterio.open(path) as src:
data = src.read(1)
# 忽略无效值（如0）
valid_mask = data != 0
total_pixels = valid_mask.sum()
area_per_pixel = 10 * 10 # 10m × 10m = 100 m²
crop_areas = {}
for crop_id in [1, 2, 3]:
crop_mask = (data == crop_id) & valid_mask
area_m2 = crop_mask.sum() * area_per_pixel
crop_areas[f"crop_{crop_id}"] = area_m2 / 1e6 # 转换为平方公里
results.append({"year": year, **crop_areas})
df = pd.DataFrame(results)
print(df)
结果示例：
| year | crop_1 | crop_2 | crop_3 |
|---|---|---|---|
| 2017 | 120.5 | 850.3 | 320.1 |
| 2018 | 125.2 | 848.7 | 325.4 |
| ... | ... | ... | ... |
该表可用于绘制趋势图，揭示作物种植结构演变规律。
可视化呈现：让数据"说话"
借助Matplotlib生成年度作物分布热力图对比：
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 4, figsize=(16, 8))
axes = axes.ravel()
for idx, year in enumerate(range(2017, 2024)):
with rasterio.open(f"data/NE_Crop_{year}.tif") as src:
data = src.read(1)
ax = axes[idx]
im = ax.imshow(data, cmap='tab20', vmin=0, vmax=3)
ax.set_title(f"{year}")
ax.axis('off')
plt.colorbar(im, ax=axes, label='作物类型')
plt.tight_layout()
plt.show()
此图直观展示十年间作物格局变迁，适用于报告展示与政策解读。
关键总结
- 所有处理均基于标准地理空间库，具备良好兼容性。
- 使用窗口读取和内存管理策略应对大文件。
- 统一坐标系是实现时间序列分析的前提。
- 面积计算结合像素计数与空间分辨率，保证精度。
- 可视化环节强化了数据解释力。
通过这套流程，原本难以驾驭的10米级高维数据，已转化为清晰、可信的农业决策支持依据。
