土地利用变化分析实战：用Python处理40年CNLUCC数据集

土地利用变化分析实战用Python处理40年CNLUCC数据集1972年至今的中国土地利用变化数据如同一部记录国土变迁的生态相册。对于区域规划师、生态研究者而言这套CNLUCC数据集的价值不亚于考古学家手中的碳14检测仪。本文将带您用Python打开这本相册从数据清洗到时空模式挖掘一步步解码40年间土地演变的秘密。1. 环境配置与数据准备工欲善其事必先利其器。处理地理空间数据需要特定的工具链组合# 推荐使用conda创建专用环境 conda create -n landuse python3.9 conda install -c conda-forge gdal geopandas matplotlib scikit-learn pip install earthpy rasterioCNLUCC数据通常以GeoTIFF格式分发每个年份对应一个文件。数据组织结构示例CNLUCC_1972.tif CNLUCC_1985.tif ... CNLUCC_2023.tif注意早期数据(1972-1985)是多年合成的结果1986年后为年度数据。分析时需注意这种时间分辨率的差异。数据加载基础代码框架import geopandas as gpd import rasterio def load_landuse(year): with rasterio.open(fCNLUCC_{year}.tif) as src: data src.read(1) meta src.meta return data, meta2. 数据清洗与标准化原始数据往往存在需要处理的常见问题边缘区域的NoData值不同年份间分类体系的微小差异投影系统的验证数据质量检查清单检查各年份数据的投影是否一致验证分类编码是否连续统计各类型面积变化曲线是否合理检查影像边缘是否有异常值投影验证代码示例from pyproj import CRS def check_projection(file_path): with rasterio.open(file_path) as src: crs CRS.from_dict(src.crs) assert crs CRS.from_proj4(projaea lat_00 lon_0105 lat_125 lat_247 x_00 y_00 ellpskrass unitsm no_defs)3. 变化检测算法实现3.1 转移矩阵分析土地类型变化的护照记录import numpy as np from sklearn.metrics import confusion_matrix def transition_matrix(data1, data2, n_classes25): # 展平数据 flat1 data1.flatten() flat2 data2.flatten() # 计算转移矩阵 mat confusion_matrix(flat1, flat2, labelsrange(n_classes)) return mat典型输出示例1972-2023年某区域类型耕地林地水域建设用地未利用地耕地65%12%3%18%2%林地8%75%2%5%10%3.2 变化热点识别使用卷积运算检测变化密集区from scipy.ndimage import uniform_filter def change_density_map(change_map, window_size5): # 二值化变化图 binary_change (change_map 0).astype(int) # 滑动窗口统计 density uniform_filter(binary_change, sizewindow_size) return density4. 时空可视化技巧4.1 动态变化图制作import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def create_animation(years): fig, ax plt.subplots(figsize(10, 8)) def update(frame): ax.clear() data, _ load_landuse(years[frame]) im ax.imshow(data, cmaptab20, vmin1, vmax25) ax.set_title(fLand Use {years[frame]}) return im ani FuncAnimation(fig, update, frameslen(years), interval500) return ani4.2 三维时空立方体展示from mpl_toolkits.mplot3d import Axes3D def plot_3d_timeseries(data_cube): fig plt.figure(figsize(12, 10)) ax fig.add_subplot(111, projection3d) # 为每类创建分层 for class_id in range(1, 26): mask (data_cube class_id) z, y, x np.where(mask) ax.scatter(x, y, z, s1, labelfClass {class_id}) ax.set_zlabel(Year) ax.legend()5. 变化模式挖掘实战5.1 驱动因子分析构建土地利用变化与社会经济指标的关联模型import statsmodels.api as sm def driver_analysis(change_rate, factors): # factors: 二维数组每列代表一个驱动因子(如GDP、人口等) X sm.add_constant(factors) model sm.OLS(change_rate, X) results model.fit() return results.summary()典型驱动因子可能包括到城市中心的距离高程和坡度道路密度人口增长率产业政策指标5.2 预测模型构建基于随机森林的变化预测from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split def train_change_model(prev_year, next_year, drivers): # 准备训练数据 X drivers.reshape(-1, drivers.shape[-1]) y (next_year ! prev_year).astype(int).flatten() # 训练测试分割 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3) # 训练模型 model RandomForestClassifier(n_estimators100) model.fit(X_train, y_train) return model, model.score(X_test, y_test)6. 专题分析案例城市扩张生态影响以长三角城市群为例展示如何量化建设用地扩张对生态系统的压力def ecological_pressure(index): 计算生态压力指数 :param index: 土地利用类型索引图 :return: 压力指数矩阵 pressure_weights { 1: 0.2, # 耕地 2: 0.1, # 林地 3: 0.3, # 草地 4: 1.0, # 建设用地 5: 0.5 # 未利用地 } pressure np.zeros_like(index, dtypefloat) for code, weight in pressure_weights.items(): pressure[index code] weight return pressure典型分析流程提取各时期建设用地分布计算缓冲区生态压力指数分析压力与生物多样性指标的相关系数识别生态敏感区变化趋势在最近的一个区域规划项目中使用这种方法发现某新城规划区将影响3个重要湿地斑块的连通性促使规划方案调整保留了生态廊道。