从Landsat到CLCD：用Python量化中国城市扩张的实战指南
站在北京国贸三期俯瞰这座城市的天际线，玻璃幕墙的摩天大楼如同森林般拔地而起。作为一名城市规划研究者，我常常思考：过去三十年间，这座城市的物理边界究竟扩张了多少？钢筋混凝土究竟吞噬了多少农田和绿地？直到接触到武汉大学黄昕团队发布的CLCD数据集，这些问题才有了精确量化的可能。本文将带你用Python从零开始，完成城市扩张分析的完整流程——从数据下载到空间统计，从趋势可视化到专业图表输出。
1. 环境配置与数据准备
工欲善其事，必先利其器。在开始分析前，我们需要搭建一个高效的Python地理数据处理环境。推荐使用conda创建专属的虚拟环境：
conda create -n urban_growth python=3.9 conda activate urban_growth conda install -c conda-forge geopandas rasterio matplotlib numpy scipy jupyterlabCLCD数据集提供了1985-2023年中国逐年30米分辨率的土地覆盖信息，其中不透水面(impervious surface)数据正是我们分析城市扩张的关键指标。数据下载可通过以下两种方式：
- 官方渠道：访问Zenodo数据仓库(doi:10.5281/zenodo.4417809)获取完整年度数据集
- 预处理版本：第三方整理的按省份分割版本(需注意数据时效性)
下载后的GeoTIFF文件建议按以下结构组织：
CLCD_data/ ├── raw/ │ ├── CLCD_v01_1990.tif │ ├── CLCD_v01_2000.tif │ └── ... └── processed/2. 数据读取与预处理
CLCD数据采用数字编码表示不同土地类型，其中不透水面对应的值为：
- 1：低密度不透水面(如居民区)
- 2：高密度不透水面(如商业区)
使用rasterio读取数据时，需要注意中国区域的数据投影通常为WGS84 UTM Zone 49N(EPSG:32649)：
import rasterio def load_clcd(year, data_dir): file_path = f"{data_dir}/raw/CLCD_v01_{year}.tif" with rasterio.open(file_path) as src: data = src.read(1) profile = src.profile return data, profile # 示例：读取2000年数据 data_2000, meta_2000 = load_clcd(2000, "CLCD_data")为提取特定城市区域，我们需要准备行政边界矢量数据。推荐使用高精度的城市行政区划shp文件：
import geopandas as gpd def load_city_boundary(city_name): # 假设已有城市边界shp文件 gdf = gpd.read_file(f"boundary/{city_name}.shp") # 统一转换为与CLCD相同的CRS return gdf.to_crs(epsg=32649) beijing = load_city_boundary("Beijing")3. 城市扩张量化分析
3.1 不透水面提取与统计
通过空间掩模运算，我们可以精确计算城市不透水面的面积变化：
import numpy as np from rasterio.mask import mask def extract_impervious(data, transform, city_geom): # 创建不透水面掩膜(值为1或2) impervious_mask = np.isin(data, [1, 2]) # 计算总像素数 total_pixels = np.sum(impervious_mask) # 转换为平方公里(30m分辨率) return total_pixels * 0.03 * 0.03 # 时间序列分析示例 years = range(1990, 2021, 5) results = [] for year in years: data, meta = load_clcd(year, "CLCD_data") area = extract_impervious(data, meta['transform'], beijing.geometry[0]) results.append((year, area))3.2 扩张趋势可视化
将统计结果转化为直观的图表：
import matplotlib.pyplot as plt # 准备数据 years = [r[0] for r in results] areas = [r[1] for r in results] # 创建图表 plt.figure(figsize=(10, 6)) plt.plot(years, areas, 'bo-', linewidth=2) plt.title('北京不透水面面积变化 (1990-2020)', fontsize=14) plt.xlabel('年份', fontsize=12) plt.ylabel('面积 (km²)', fontsize=12) plt.grid(True, alpha=0.3) plt.savefig('beijing_growth.png', dpi=300, bbox_inches='tight')典型城市扩张分析可关注以下指标：
| 指标名称 | 计算公式 | 生态意义 |
|---|---|---|
| 年扩张速率 | (A₂-A₁)/(A₁×Δt) | 城市发展强度 |
| 扩张弹性系数 | 城市面积增速/人口增速 | 土地集约利用程度 |
| 景观破碎度指数 | 边界长度/面积 | 城市形态规整度 |
4. 空间格局演变分析
4.1 扩张方向识别
通过扇形分析法可以量化城市扩张的方向偏好：
def sector_analysis(city_center, years, radius=50000, n_sectors=8): angles = np.linspace(0, 360, n_sectors+1) growth_by_sector = {i: [] for i in range(n_sectors)} for year in years: data, _ = load_clcd(year, "CLCD_data") # 实现扇形区域统计逻辑... return growth_by_sector4.2 空间热点检测
使用Getis-Ord Gi*统计量识别扩张热点区域：
from pysal.explore import esda from pysal.lib import weights def detect_hotspots(current_year, previous_year): # 计算变化网格 diff = current_year - previous_year # 创建空间权重矩阵 w = weights.DistanceBand.from_array(diff, threshold=3000) # 计算Gi*统计量 gi = esda.G_Local(diff, w) return gi.z_sim5. 高级分析与应用
5.1 驱动因子分析
结合社会经济数据探究城市扩张驱动力：
import pandas as pd from sklearn.ensemble import RandomForestRegressor def driving_force_analysis(city_name): # 准备数据 growth_df = pd.read_csv(f"{city_name}_growth.csv") economy_df = pd.read_csv(f"{city_name}_gdp.csv") # 数据合并与特征工程 merged = pd.merge(growth_df, economy_df, on='year') X = merged[['gdp', 'population', 'fdi']] y = merged['growth_rate'] # 训练模型 model = RandomForestRegressor() model.fit(X, y) return model.feature_importances_5.2 未来情景模拟
基于历史趋势预测城市发展：
from statsmodels.tsa.arima.model import ARIMA def predict_growth(history_years, history_areas, pred_years=5): model = ARIMA(history_areas, order=(1,1,1)) model_fit = model.fit() forecast = model_fit.forecast(steps=pred_years) return forecast在完成上海浦东新区的分析项目时，我发现2005-2010年间的不透水面增长出现了异常波动。核查原始影像后发现，这一时期该区域正在进行大规模填海造地工程，导致CLCD分类出现短暂偏差。这种实地验证的过程让我深刻认识到：任何自动化分类结果都需要结合地面真实情况和历史背景进行解读。
