中国土地利用30年变迁可视化指南：基于CLCD数据的Python分析全流程

中国土地利用30年变迁可视化指南：基于CLCD数据的Python分析全流程

过去三十年，中国经历了人类历史上规模最大、速度最快的城市化进程。从沿海经济特区的崛起到西部大开发的推进，每一寸土地的变迁都承载着经济发展的密码。对于数据爱好者和城市规划研究者而言，如何将这些变化转化为直观的可视化呈现，不仅是一项技术挑战，更是一次发现中国发展脉络的探索之旅。

本文将带您从零开始，使用Python技术栈处理CLCD（中国土地利用覆盖数据集）的TIFF格式数据，通过完整的工作流展示从数据获取到时空分析，再到动态可视化呈现的全过程。不同于简单的数据介绍，我们将重点关注如何挖掘数据背后的故事，让冰冷的数字讲述生动的变迁历程。

1. 环境准备与数据获取

1.1 Python环境配置

处理地理空间数据需要特定的Python库支持。推荐使用Anaconda创建独立环境，避免与其他项目产生依赖冲突：

conda create -n landuse python=3.9 conda activate landuse pip install rasterio geopandas matplotlib numpy pandas xarray cartopy ipyleaflet

关键库说明：

• rasterio：专业的栅格数据处理库
• geopandas：地理空间数据分析的利器
• cartopy：专业的地图投影与可视化工具
• ipyleaflet：交互式地图展示

提示：如果安装cartopy遇到问题，可先通过conda安装基础依赖：conda install -c conda-forge cartopy

1.2 CLCD数据下载与初步检查

CLCD数据集可从武汉大学杨杰教授团队提供的公开渠道获取。数据集按年份组织，每个年份包含中国全境的30米分辨率土地利用分类TIFF文件。

下载完成后，建议先进行数据完整性检查：

import rasterio def check_tiff(filepath): try: with rasterio.open(filepath) as src: print(f"文件有效，基本信息：") print(f"尺寸：{src.width}x{src.height}") print(f"波段数：{src.count}") print(f"坐标参考系统(CRS)：{src.crs}") print(f"边界范围：{src.bounds}") return True except Exception as e: print(f"文件损坏或格式错误：{str(e)}") return False # 示例检查1990年数据 check_tiff("CLCD_1990.tif")

常见问题排查表：

2. 数据处理与质量评估

2.1 数据预处理流程

原始TIFF数据需要经过一系列预处理才能用于分析。以下是标准化的预处理流程：

1. 坐标系统一化：确保所有年份数据使用相同的CRS（WGS84）
2. 无效值处理：识别并处理Nodata值
3. 边界对齐：保证不同年份数据像素对齐
4. 内存优化：分块处理大数据量

import numpy as np from rasterio.windows import Window def preprocess_landuse(input_path, output_path): with rasterio.open(input_path) as src: # 创建处理参数配置 profile = src.profile profile.update( compress='lzw', nodata=0 ) # 分块处理大文件 block_size = 2048 height, width = src.shape with rasterio.open(output_path, 'w', **profile) as dst: for i in range(0, height, block_size): for j in range(0, width, block_size): window = Window( j, i, min(block_size, width - j), min(block_size, height - i) ) data = src.read(window=window) # 简单清洗：将异常值转为0 data[data > 8] = 0 dst.write(data, window=window)

2.2 数据质量评估指标

CLCD数据集虽然质量较高，但仍需进行质量评估。我们可以从三个维度进行验证：

1. 时间一致性检查

import xarray as xr def check_temporal_consistency(files): # 创建时间序列数据集 datasets = [] years = [] for f in sorted(files): year = int(f.split('_')[1].split('.')[0]) with rasterio.open(f) as src: data = src.read(1) datasets.append(data) years.append(year) # 计算年际变化率 changes = [] for i in range(1, len(datasets)): change = np.sum(datasets[i] != datasets[i-1]) / datasets[i].size changes.append(change) return years[1:], changes

2. 分类准确性验证

通过与公开的抽样验证点对比，计算混淆矩阵：

3. 空间自相关分析

使用Moran's I指数检查空间分布模式：

from libpysal.weights import lat2W from esda.moran import Moran def calculate_spatial_autocorrelation(data): # 创建空间权重矩阵 w = lat2W(data.shape[0], data.shape[1]) # 计算Moran's I moran = Moran(data.ravel(), w) return moran.I, moran.p_sim

3. 变化检测与时空分析

3.1 土地利用转移矩阵

土地利用变化的核心是类型间的转换关系。转移矩阵能清晰展示这种动态：

def create_transition_matrix(from_year, to_year): with rasterio.open(f'CLCD_{from_year}.tif') as src1, \ rasterio.open(f'CLCD_{to_year}.tif') as src2: data1 = src1.read(1) data2 = src2.read(1) # 创建9x9转移矩阵 matrix = np.zeros((9, 9), dtype=np.int64) for i in range(9): mask = (data1 == i+1) for j in range(9): matrix[i, j] = np.sum((data2 == j+1) & mask) return matrix

示例输出（1990-2020年部分区域）：

3.2 热点区域识别

使用Getis-Ord Gi*统计量识别变化热点：

import pysal as ps from esda.getisord import G_Local def detect_hotspots(change_layer): # 标准化变化强度 zscore = (change_layer - np.mean(change_layer)) / np.std(change_layer) # 创建空间权重 w = ps.weights.lat2W(change_layer.shape[0], change_layer.shape[1]) # 计算Gi* gi = G_Local(zscore.flatten(), w) # 重塑结果 hotspots = gi.Zs.reshape(change_layer.shape) return hotspots

热点分级标准：

4. 高级可视化技术

4.1 动态时空演变展示

使用ipyleaflet创建交互式时间轴地图：

from ipyleaflet import Map, basemaps, GeoJSON import ipywidgets as widgets def create_time_slider_map(bbox): # 初始化地图 m = Map( center=((bbox[1]+bbox[3])/2, (bbox[0]+bbox[2])/2), zoom=6, basemap=basemaps.Esri.WorldImagery ) # 创建时间滑块 year_slider = widgets.IntSlider( value=1990, min=1990, max=2020, step=5, description='年份:' ) # 颜色映射 cmap = { 1: '#FAFAD2', # 耕地 2: '#228B22', # 林地 3: '#9ACD32', # 灌木 4: '#7CFC00', # 草地 5: '#1E90FF', # 水域 6: '#FFFFFF', # 冰雪 7: '#CD853F', # 裸地 8: '#FF0000', # 建设用地 9: '#9370DB' # 湿地 } # 交互更新函数 def update_map(year): with rasterio.open(f'CLCD_{year}.tif') as src: data = src.read(1) # 转换为GeoJSON格式 shapes = rasterio.features.shapes(data, transform=src.transform) features = [] for geom, val in shapes: if val > 0: features.append({ 'type': 'Feature', 'geometry': geom, 'properties': {'value': int(val)} }) geojson = { 'type': 'FeatureCollection', 'features': features } # 更新地图 if hasattr(m, 'landuse_layer'): m.remove_layer(m.landuse_layer) m.landuse_layer = GeoJSON( data=geojson, style={ 'color': 'black', 'opacity': 0.2, 'fillOpacity': 0.7, 'fillColor': { 'property': 'value', 'type': 'categorical', 'values': list(cmap.keys()), 'colors': list(cmap.values()) } } ) m.add_layer(m.landuse_layer) # 绑定交互 widgets.interactive(update_map, year=year_slider) # 组合控件 controls = widgets.HBox([year_slider]) return widgets.VBox([controls, m])

4.2 三维变化地形图

使用matplotlib创建三维地表变化视图：

from mpl_toolkits.mplot3d import Axes3D def plot_3d_change(start_year, end_year, region): fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') # 加载数据 with rasterio.open(f'CLCD_{start_year}.tif') as src1, \ rasterio.open(f'CLCD_{end_year}.tif') as src2: # 裁剪感兴趣区域 window = rasterio.windows.from_bounds(*region, transform=src1.transform) data1 = src1.read(1, window=window) data2 = src2.read(1, window=window) # 计算变化矩阵 change = np.where(data1 != data2, 1, 0) # 创建网格 rows, cols = np.indices(change.shape) x = cols.flatten() y = rows.flatten() z = np.zeros_like(x) dx = dy = np.ones_like(x) dz = change.flatten() * 50 # 放大变化高度 # 绘制柱状图 ax.bar3d(x, y, z, dx, dy, dz, color='r', alpha=0.3) # 设置视角 ax.view_init(elev=45, azim=45) ax.set_title(f'土地利用变化三维视图 {start_year}-{end_year}') ax.set_xlabel('经度方向') ax.set_ylabel('纬度方向') ax.set_zlabel('变化强度') plt.tight_layout() return fig

在实际项目中，我发现将三维可视化与时间滑块结合，能更直观地展示城市扩张的时空模式。特别是对于长三角、珠三角等快速城市化区域，这种可视化方式能清晰呈现"摊大饼"式的扩张特征。