城市扩张可视化用Python解码30年不透水层变迁故事当一张张卫星影像从高空俯瞰大地那些灰白色的斑块如同城市的年轮记录着人类文明扩张的足迹。这些被称为“不透水层”的区域——建筑、道路、广场等人工硬化地表是城市化进程最直观的物理印记。对于城市规划研究者、环境科学家乃至政策制定者而言如何从海量的遥感数据中提取出有温度、有深度的城市发展叙事是一个既专业又充满挑战的课题。传统的城市规划分析往往依赖于统计年鉴和宏观数据缺乏时空上的精细刻画。而如今像GAIAGlobal Artificial Impervious Area这样的全球逐年不透水层数据集为我们打开了一扇窗。它提供了自1985年以来全球范围内30米分辨率的逐年城市地表覆盖信息。每一个像素值都像是一个时间胶囊封装了城市在某一年份出现或持续存在的“记忆”。然而这些数据以冷冰冰的GeoTIFF格式存在如何用Python将它们转化为动态的增长热力图并与人口、经济等数据交叉融合讲出一个既有科研价值又有传播力的时空故事这正是本文要探讨的核心。我们将超越简单的数据读取与静态出图深入一套完整的分析流程从多时相TIFF数据的批量处理与异常值清洗到城市扩张热点的时空检测与可视化表达再到结合外部数据进行驱动因子交叉验证。整个过程我们将使用Python生态中的rasterio、xarray、geopandas、matplotlib以及plotly等库构建一个可复现、可扩展的分析框架。无论你是刚接触空间数据分析的研究生还是希望将研究成果可视化的资深学者都能从中获得实用的代码范例和创新的分析视角。1. 数据基石理解与准备全球不透水层时序数据在开始编码之前我们必须先理解手中“武器”的特性。GAIA数据集并非一个简单的全球单文件而是采用了一种分幅存储的策略。它将全球纬度范围-60°到80°划分为5°×5°的渔网格网每个格网保存为一个独立的TIFF文件。这种设计兼顾了数据管理的便利性与局部区域分析的灵活性。文件命名规则是数据访问的钥匙。每个文件遵循GAIA_1985_2024_{经度}_{纬度}.tif的格式其中的经度和纬度指的是该5°×5°格网左上角的坐标。例如GAIA_1985_2024_110_25.tif表示其覆盖范围是东经110°至115°北纬25°至20°注意纬度递减。此外每个图幅边缘包含了30米的缓冲区确保了相邻图幅在镶嵌时能够无缝衔接。注意在下载数据时你需要根据研究区域的范围计算出覆盖该区域所需的所有格网文件。例如研究整个长江三角洲可能需要下载多个相邻的格网文件并进行镶嵌处理。数据最核心的信息蕴藏在像元值中。GAIA采用了一种巧妙的编码方式像元值范围是0-40。0: 代表该像元在1985-2024年间从未被识别为城市区域。1-40: 代表城市区域首次被识别出的年份。1对应2024年新增2对应2023年新增以此类推40对应1985年及之前就已存在的城市区域。这意味着一张TIFF图像实际上压缩了40年的城市变迁历史。我们的任务就是将其“解压缩”并展开成时间序列。为了高效处理这些可能多达数十甚至上百个的TIFF文件我们首先需要建立一个清晰的项目目录结构和数据加载管道。以下是一个推荐的项目初始化代码它定义了数据路径并创建了必要的文件夹。import os from pathlib import Path # 定义项目根目录和子目录 PROJECT_ROOT Path(/path/to/your/project) # 请替换为你的实际路径 DATA_RAW PROJECT_ROOT / data/raw_gaia # 存放原始下载的GAIA分幅TIFF DATA_PROCESSED PROJECT_ROOT / data/processed # 存放处理后的数据如镶嵌、重采样结果 DATA_EXTERNAL PROJECT_ROOT / data/external # 存放人口、GDP等外部数据 RESULTS PROJECT_ROOT / results # 存放最终图表、报告 SCRIPTS PROJECT_ROOT / scripts # 存放分析脚本 # 创建目录如果不存在 for dir_path in [DATA_RAW, DATA_PROCESSED, DATA_EXTERNAL, RESULTS, SCRIPTS]: dir_path.mkdir(parentsTrue, exist_okTrue) print(项目目录结构初始化完成。)2. 核心引擎用Python高效读取与预处理时序栅格有了数据下一步就是将其读入Python环境并进行初步的“体检”与清洗。我们将主要使用rasterio和xarray库。rasterio提供了丰富的地理空间数据读写接口而xarray特别擅长处理多维数组如时间序列的栅格数据并集成了dask以实现并行计算非常适合处理大数据量的遥感影像。2.1 单时相数据读取与坐标信息提取首先我们来看如何读取单个年份实际上是包含所有年份信息的一个文件的数据并正确提取其地理坐标信息。与原文仅展示绘图不同我们更关注数据的元信息这是后续所有时空分析的基础。import rasterio import numpy as np import xarray as xr # 示例读取一个GAIA分幅文件 tiff_path DATA_RAW / GAIA_1985_2024_110_25.tif with rasterio.open(tiff_path) as src: # 读取第一个波段的数据GAIA数据通常只有一个波段 data src.read(1) # 返回一个二维numpy数组 # 获取关键的元数据 meta src.meta.copy() transform src.transform # 仿射变换矩阵用于计算地理坐标 crs src.crs # 坐标参考系统GAIA通常是WGS84 (EPSG:4326) bounds src.bounds # 图像的地理边界 (left, bottom, right, top) height, width src.height, src.width print(f数据形状: {data.shape}) print(f数据类型: {data.dtype}) print(f坐标参考系: {crs}) print(f地理范围: {bounds}) print(f仿射变换参数: {transform})这段代码不仅读出了像元值矩阵更重要的是获取了将像素坐标转换为真实世界坐标经纬度所需的全部信息。transform是一个包含6个参数的元组(c, a, b, f, d, e)其中c和f是左上角像素的X经度和Y纬度坐标a和e是像素在X和Y方向上的大小。2.2 批量读取与时间维度构建真正的挑战在于处理多年份、多图幅的数据。我们需要将每个分幅文件中的“年份编码”数据转换成一个清晰的时间序列数据立方体Data Cube。思路是遍历所有相关分幅文件根据像元值1-40反向推算出每个年份的城市分布二值图0或1然后将所有年份的数据在时间维度上堆叠起来。下面的函数演示了如何将单个GAIA文件解码为一个xarray.DataArray其中包含一个名为year的维度。def decode_gaia_to_timeseries(tiff_path, start_year1985, end_year2024): 将GAIA TIFF文件解码为时间序列的xarray DataArray。 参数: tiff_path: GAIA TIFF文件路径。 start_year, end_year: 数据的时间范围。 返回: xarray.DataArray维度为 (year, y, x)。 with rasterio.open(tiff_path) as src: data src.read(1) transform src.transform crs src.crs height, width data.shape # 生成年份列表 years np.arange(start_year, end_year 1) # 初始化一个三维数组来存储每年的二值图 yearly_data np.zeros((len(years), height, width), dtypenp.uint8) # 核心解码逻辑根据像元值v将其标记为对应年份及之前所有年份的城市区域 for v in range(1, 41): # v从1到40 target_year end_year - v 1 # 计算该值对应的年份 year_index np.where(years target_year)[0][0] # 找出所有像元值 v 的像素这些像素在target_year及之后都是城市 # 但为了得到每年独立的图我们标记该年份“存在”城市 mask (data v) (data 40) # 更精确的做法标记该年份“新增”的城市。这里简化处理标记为“存在” yearly_data[year_index:, mask] 1 # 从该年份开始往后都标记为1 # 创建xarray DataArray lon np.linspace(transform.c, transform.c transform.a * width, width) lat np.linspace(transform.f, transform.f transform.e * height, height) da xr.DataArray( yearly_data, dims[year, lat, lon], coords{ year: years, lat: lat, lon: lon }, attrs{ crs: crs, transform: transform, description: fDecoded urban impervious area from {tiff_path.name} } ) return da # 使用示例 decoded_da decode_gaia_to_timeseries(tiff_path) print(decoded_da)这个函数是理解GAIA数据时间维度的关键。它输出的DataArray是一个三维数据立方体我们可以轻松地对其进行时间切片如decoded_da.sel(year2010)获取2010年的城市分布图或者进行时间序列的统计分析。2.3 多图幅镶嵌与异常值处理研究区域往往跨越多个5°×5°的格网。我们需要将这些分幅数据“缝合”起来形成一个完整的研究区覆盖。这个过程称为镶嵌Mosaicking。同时原始数据中可能存在噪声或异常值如因云层、传感器故障导致的错误分类需要进行清洗。多图幅镶嵌可以使用rasterio.merge模块完成。我们需要先找到所有覆盖研究区的文件然后进行合并。from rasterio.merge import merge import glob # 假设研究区中心在(经度, 纬度)我们找出周围9个格网 center_lon, center_lat 112.5, 23.5 # 计算格网左上角坐标5度的整数倍 grid_lon (center_lon // 5) * 5 grid_lat ((center_lat // 5) 1) * 5 # 纬度向上取整 tiff_files [] for lon_offset in [-5, 0, 5]: for lat_offset in [-5, 0, 5]: target_lon grid_lon lon_offset target_lat grid_lat lat_offset pattern fGAIA_1985_2024_{int(target_lon)}_{int(target_lat)}.tif file_path DATA_RAW / pattern if file_path.exists(): tiff_files.append(file_path) print(f找到 {len(tiff_files)} 个需要镶嵌的文件。) # 执行镶嵌 src_files_to_mosaic [rasterio.open(fp) for fp in tiff_files] mosaic, out_transform merge(src_files_to_mosaic) # 关闭所有文件 for src in src_files_to_mosaic: src.close() # 保存镶嵌结果 output_mosaic_path DATA_PROCESSED / study_area_mosaic.tif with rasterio.open( output_mosaic_path, w, driverGTiff, heightmosaic.shape[1], widthmosaic.shape[2], count1, dtypemosaic.dtype, crssrc_files_to_mosaic[0].crs, transformout_transform, ) as dst: dst.write(mosaic)异常值处理通常结合数据本身的特性。对于GAIA数据像元值理论上应在0-40之间。我们可以设置一个简单的阈值进行过滤。更高级的方法可能涉及邻域分析或与更高精度的参考数据对比但这超出了基础处理的范畴。# 简单的异常值过滤将不在0-40范围内的值设为无数据例如-9999 decoded_da_clean decoded_da.where((decoded_da 0) (decoded_da 40), other-9999) # 或者如果我们只关心城市/非城市可以将异常值归为“非城市”(0) decoded_da_binary xr.where((decoded_da 1) (decoded_da 40), 1, 0)3. 时空叙事从静态出图到动态分析与热点探测将数据处理妥当后就进入了最富创造性的环节——可视化与叙事。我们的目标不仅是画出某一年份的地图更是要揭示城市扩张的动态过程、空间模式及其背后的驱动因素。3.1 制作城市扩张动态热力图静态地图难以表现“变迁”。我们可以计算每个像元在特定时间段内例如每5年或10年从非城市变为城市的频率或概率生成一张“扩张强度”热力图。这比单纯展示两个时间点的差异图包含了更多信息。import matplotlib.pyplot as plt import matplotlib.colors as mcolors import cartopy.crs as ccrs import cartopy.feature as cfeature # 假设我们已经有了解码后的时间序列数据立方体 urban_cube (year, lat, lon) # 计算1985-2000年间的扩张总量新增城市像元数 urban_start urban_cube.sel(year1985) urban_end urban_cube.sel(year2000) expansion_85_00 urban_end - urban_start # 结果为1扩张、0不变、-1萎缩理论上GAIA数据中极少 # 或者计算扩张面积需要知道每个像元的实际面积在WGS84经纬度下需用球面面积公式近似 # 这里我们简单地将扩张像元数作为强度 expansion_intensity expansion_85_00.where(expansion_85_00 1, 0) # 只保留扩张区域 # 绘制扩张热力图 fig plt.figure(figsize(14, 10)) ax plt.axes(projectionccrs.PlateCarree()) ax.set_extent([110, 118, 20, 26]) # 设定绘图范围例如珠三角地区 # 添加地理背景 ax.add_feature(cfeature.COASTLINE.with_scale(50m), linewidth0.8) ax.add_feature(cfeature.BORDERS.with_scale(50m), linestyle:, linewidth0.5) ax.add_feature(cfeature.RIVERS.with_scale(50m), linewidth0.5, alpha0.5) # 绘制扩张强度 im ax.imshow(expansion_intensity.squeeze(), # 去掉可能的单维度 extent[urban_cube.lon.min(), urban_cube.lon.max(), urban_cube.lat.min(), urban_cube.lat.max()], originupper, cmapYlOrRd, # 黄-橙-红色系适合表示强度 transformccrs.PlateCarree(), interpolationnearest) # 添加颜色条 cbar plt.colorbar(im, axax, orientationvertical, pad0.05, shrink0.7) cbar.set_label(Urban Expansion Intensity (1985-2000), fontsize12) ax.set_title(Urban Expansion Hotspots (1985-2000), fontsize16, fontweightbold) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) plt.tight_layout() plt.savefig(RESULTS / expansion_hotspots_85_00.png, dpi300, bbox_inchestight) plt.show()3.2 时间序列折线图与统计指标除了空间模式时间趋势同样重要。我们可以选取整个研究区域或者某个特定的行政区如一个城市计算其每年的不透水层总面积绘制时间序列曲线。# 计算研究区每年城市像元的总数假设每个像元代表城市 # 注意这并非实际面积要计算面积需要投影转换。这里用像元数作为相对指标。 annual_urban_pixels urban_cube.sum(dim[lat, lon]) # 对空间维度求和 # 转换为DataFrame便于绘图和分析 import pandas as pd urban_ts_df annual_urban_pixels.to_dataframe(nameurban_pixel_count).reset_index() plt.figure(figsize(12, 6)) plt.plot(urban_ts_df[year], urban_ts_df[urban_pixel_count], markero, linewidth2, markersize4) plt.fill_between(urban_ts_df[year], urban_ts_df[urban_pixel_count], alpha0.3) plt.xlabel(Year, fontsize12) plt.ylabel(Urban Pixel Count (Relative Area), fontsize12) plt.title(Temporal Trend of Impervious Surface Area (1985-2024), fontsize14, fontweightbold) plt.grid(True, linestyle--, alpha0.7) plt.xticks(urban_ts_df[year][::5]) # 每5年显示一个刻度 plt.tight_layout() plt.savefig(RESULTS / urban_area_trend.png, dpi300) plt.show()为了更深入地分析我们可以计算一些关键的扩张指标指标公式示例描述年均扩张速率(Area_end - Area_start) / (Year_end - Year_start)衡量扩张的绝对速度年均扩张强度指数(Area_end - Area_start) / (Area_start * (Year_end - Year_start)) * 100%衡量相对于基期面积的扩张强度扩张的空间紧凑度(Perimeter^2) / (4π * Area)形状指数值越大越不规则扩张可能越分散重心迁移轨迹计算每年城市像元的平均经纬度可视化城市扩张的主要方向这些指标的计算可以通过xarray的空间运算和geopandas的几何操作相结合来实现。3.3 交互式可视化与故事讲述静态图表适合论文报告而交互式图表则能极大地增强数据的探索性和故事的感染力。plotly或folium库可以帮助我们创建可缩放、可悬停查看详细信息的地图。下面是一个使用plotly创建交互式时间滑动条的示例可以动态展示城市逐年扩张的过程。import plotly.express as px import plotly.graph_objects as go from plotly.subplots import make_subplots # 为了性能可能需要对数据进行重采样降低空间分辨率和年份抽样 urban_cube_sampled urban_cube.sel(yearslice(1985, 2024, 5)) # 每5年取一帧 # 将xarray DataArray转换为适合plotly的格式这里需要一些数据重塑操作 # 假设我们已将数据转换为一个DataFrame包含列lon, lat, year, is_urban # df_urban_long ... (数据转换过程略) # 创建动画地图 fig px.scatter_mapbox(df_urban_long[df_urban_long[is_urban]1], latlat, lonlon, animation_frameyear, # 时间轴 color_discrete_sequence[red], # 城市区域用红色表示 zoom8, height600, titleUrban Expansion Animation (1985-2024, 5-year interval)) fig.update_layout(mapbox_styleopen-street-map) # 使用OpenStreetMap底图 fig.update_layout(margin{r:0,t:30,l:0,b:0}) # 将动画保存为HTML文件可在浏览器中交互查看 fig.write_html(RESULTS / urban_expansion_animation.html)将这样的交互式地图嵌入到在线报告或演示文稿中观众可以自己拖动时间滑块亲眼目睹城市如何像生命体一样生长、蔓延这种体验远比看一系列静态图或听描述要震撼得多。4. 交叉验证与深度洞察连接空间扩张与社会经济脉搏城市扩张不是孤立的地理现象其背后是人口增长、经济发展、政策导向等多重力量的驱动。将不透水层数据与外部数据集进行空间叠加和统计分析能让我们的故事从“发生了什么”深入到“为什么发生”。4.1 与人口/经济数据的空间关联分析我们可以获取研究区域历年的人口密度栅格数据或行政区划的统计年鉴数据。通过空间统计计算每个行政区内部不透水层面积的增长与人口、GDP增长之间的相关性。操作步骤概览数据准备获取行政区划矢量边界如.shp文件和对应的社会经济数据表。空间连接使用geopandas将不透水层栅格数据按行政区进行分区统计Zonal Statistics。例如计算每个区县每年的城市像元总数。数据融合将分区统计结果与 socioeconomic 数据表通过行政区代码进行合并。相关性分析计算不透水层面积增长率与人口增长率、GDP增长率之间的皮尔逊相关系数并绘制散点图与趋势线。import geopandas as gpd import pandas as pd from rasterstats import zonal_stats # 1. 加载行政区划矢量 gdf_districts gpd.read_file(DATA_EXTERNAL / study_area_districts.shp) # 确保矢量与栅格数据在同一坐标系下WGS84 gdf_districts gdf_districts.to_crs(epsg4326) # 2. 对每一年份的不透水层二值图进行分区统计 stats_by_year [] for year in urban_cube.year.values: # 获取该年份的栅格数据需保存为临时TIFF或使用内存文件 year_data urban_cube.sel(yearyear).values # 使用rasterstats计算每个多边形内的像元统计值如总和、均值 stats zonal_stats(gdf_districts, year_data, affinetransform, stats[sum, mean]) stats_df pd.DataFrame(stats) stats_df[year] year stats_by_year.append(stats_df) # 合并所有年份的统计结果 all_stats_df pd.concat(stats_by_year, ignore_indexTrue) # 与行政区属性表合并 gdf_districts_with_stats gdf_districts.merge(all_stats_df, left_indexTrue, right_indexTrue) # 3. 加载社会经济数据假设有CSV文件包含‘district_id’, ‘year’, ‘population’, ‘gdp’ socio_df pd.read_csv(DATA_EXTERNAL / socioeconomic_data.csv) # 合并 analysis_df pd.merge(gdf_districts_with_stats, socio_df, on[district_id, year]) # 4. 计算增长率并分析相关性 # 假设我们分析2000-2020年 df_2000 analysis_df[analysis_df[year]2000].set_index(district_id) df_2020 analysis_df[analysis_df[year]2020].set_index(district_id) df_growth pd.DataFrame() df_growth[urban_area_growth] (df_2020[sum] - df_2000[sum]) / df_2000[sum] df_growth[population_growth] (df_2020[population] - df_2000[population]) / df_2000[population] df_growth[gdp_growth] (df_2020[gdp] - df_2000[gdp]) / df_2000[gdp] corr_urban_pop df_growth[urban_area_growth].corr(df_growth[population_growth]) corr_urban_gdp df_growth[urban_area_growth].corr(df_growth[gdp_growth]) print(f城市面积增长与人口增长相关系数: {corr_urban_pop:.3f}) print(f城市面积增长与GDP增长相关系数: {corr_urban_gdp:.3f})4.2 构建分析报告与叙事框架最后将所有的分析结果——动态地图、趋势图表、相关性指标——整合到一个连贯的叙事中。你可以使用Jupyter Notebook或Quarto、R Markdown通过reticulate调用Python来创建可重复生成的分析报告。报告的结构可以这样组织引言阐述研究背景、意义及使用的数据GAIA。研究方法简要说明数据处理流程解码、镶嵌、异常处理和分析方法扩张强度计算、热点探测、相关性分析。结果与分析时空格局展示研究区城市扩张的时空动态图指出扩张最快的时期和区域。扩张模式通过紧凑度、重心迁移等指标分析扩张是紧凑填充式还是蔓延式。驱动因素展示与社会经济数据的相关性分析结果讨论人口、经济与城市扩张的关系。典型案例选取一个扩张热点区域如某个新城结合卫星影像或地图进行深入解读。讨论与结论总结主要发现指出研究的局限性如GAIA数据本身的精度、社会经济数据的时空分辨率匹配问题并提出未来展望。在整个流程中我习惯将关键的参数如研究区范围、时间区间、分析指标放在配置文件如config.yaml或脚本开头这样在更换研究区域或调整分析时段时只需修改少数几个参数整个分析管道就能自动运行极大地提升了研究的可重复性和效率。处理大规模栅格数据时内存管理是个大问题xarray与dask的集成能有效进行分块处理和延迟计算对于全国甚至全球尺度的分析至关重要。