从Landsat到CLCD：手把手教你用Python分析中国城市扩张（附完整代码）

从Landsat到CLCD用Python量化中国城市扩张的实战指南站在北京国贸三期俯瞰这座城市的天际线玻璃幕墙的摩天大楼如同森林般拔地而起。作为一名城市规划研究者我常常思考过去三十年间这座城市的物理边界究竟扩张了多少钢筋混凝土究竟吞噬了多少农田和绿地直到接触到武汉大学黄昕团队发布的CLCD数据集这些问题才有了精确量化的可能。本文将带你用Python从零开始完成城市扩张分析的完整流程——从数据下载到空间统计从趋势可视化到专业图表输出。1. 环境配置与数据准备工欲善其事必先利其器。在开始分析前我们需要搭建一个高效的Python地理数据处理环境。推荐使用conda创建专属的虚拟环境conda create -n urban_growth python3.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(fboundary/{city_name}.shp) # 统一转换为与CLCD相同的CRS return gdf.to_crs(epsg32649) 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-, linewidth2) plt.title(北京不透水面面积变化 (1990-2020), fontsize14) plt.xlabel(年份, fontsize12) plt.ylabel(面积 (km²), fontsize12) plt.grid(True, alpha0.3) plt.savefig(beijing_growth.png, dpi300, bbox_inchestight)典型城市扩张分析可关注以下指标指标名称计算公式生态意义年扩张速率(A₂-A₁)/(A₁×Δt)城市发展强度扩张弹性系数城市面积增速/人口增速土地集约利用程度景观破碎度指数边界长度/面积城市形态规整度4. 空间格局演变分析4.1 扩张方向识别通过扇形分析法可以量化城市扩张的方向偏好def sector_analysis(city_center, years, radius50000, n_sectors8): angles np.linspace(0, 360, n_sectors1) 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, threshold3000) # 计算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, onyear) 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_years5): model ARIMA(history_areas, order(1,1,1)) model_fit model.fit() forecast model_fit.forecast(stepspred_years) return forecast在完成上海浦东新区的分析项目时我发现2005-2010年间的不透水面增长出现了异常波动。核查原始影像后发现这一时期该区域正在进行大规模填海造地工程导致CLCD分类出现短暂偏差。这种实地验证的过程让我深刻认识到任何自动化分类结果都需要结合地面真实情况和历史背景进行解读。