用JRC全球地表水数据集，5分钟搞定你所在城市30年水域变迁分析（附Python代码）
用JRC全球地表水数据集分析城市水域30年变迁：从数据下载到动态可视化全流程
每次路过家乡那条小时候常去的小河，总会好奇它这些年到底变宽了还是变窄了。作为地理数据爱好者，终于找到了一个专业又简单的方法——利用欧洲联合研究中心(JRC)发布的全球地表水数据集，配合Python生态中的地理处理工具，任何人都能亲手绘制出自己城市的水域变迁图。本文将带你完整走通从数据获取、区域裁剪到动态可视化的全流程，用代码重现记忆中的那片水域。
1. 认识JRC地表水数据集
JRC Global Surface Water Mapping Layers是目前最权威的全球地表水变化记录，包含1984至2019年间Landsat系列卫星拍摄的超过300万景影像。与普通卫星图不同，这套数据已经过专业处理，直接告诉我们每个30m×30m的网格是否属于水体。
数据集包含7个核心图层：
- Occurrence：每个位置在35年间被检测为水体的频率（0-100%）
- Change：前后两个时期（1984-1999 vs 2000-2019）的水体变化程度
- Extent：35年间曾出现过水体的所有位置
- Recurrence：水体每年出现的频率
- Transitions：水体状态转变（永久性/季节性/非水体间的转换）
- Seasonality：水体存在的月份数
# 数据集基本信息 import pandas as pd data_layers = { "图层名称": ["Occurrence", "Change", "Extent", "Recurrence", "Transitions", "Seasonality"], "时间范围": ["1984-2019", "1984-2019", "1984-2019", "1984-2019", "1984-2019", "2019"], "数值含义": ["水体存在概率%", "变化程度指标", "二元水体掩膜", "年际重现频率%", "状态转换编码", "存在月份数"] } pd.DataFrame(data_layers)2. 数据获取与预处理
2.1 下载指定区域数据
访问Google Earth Engine平台，注册开发者账号后，可以直接通过Python API调用数据集。为避免下载全球数据的巨大开销，我们先确定研究区域的边界坐标。
# 安装Earth Engine Python API !pip install earthengine-api # 认证并初始化 import ee ee.Authenticate() ee.Initialize() # 定义杭州市边界（示例） hangzhou = ee.FeatureCollection("users/your_account/hangzhou_boundary") # 获取2015年水体数据示例 water_2015 = ee.ImageCollection('JRC/GSW1_2/MonthlyHistory') \ .filter(ee.Filter.calendarRange(2015,2015,'year')) \ .mosaic() \ .clip(hangzhou)2.2 本地环境配置
处理地理栅格数据需要特定Python库：
# 推荐使用conda创建环境 conda create -n water_analysis python=3.8 conda activate water_analysis conda install -c conda-forge geopandas rasterio matplotlib contextily3. 城市水域变迁分析实战
3.1 行政区划裁剪技巧
以杭州市为例，我们需要从全球数据中精确提取市辖范围。建议使用GADM的行政区划数据：
import geopandas as gpd # 加载行政区划矢量边界 city_boundary = gpd.read_file("hangzhou_districts.shp") # 坐标系统一转换（关键步骤！） city_boundary = city_boundary.to_crs("EPSG:4326") # WGS84坐标系 # 可视化确认边界 ax = city_boundary.plot(edgecolor='red', facecolor='none') ctx.add_basemap(ax, crs=city_boundary.crs, source=ctx.providers.Stamen.Terrain)3.2 年度水域面积计算
通过遍历每年数据，统计水域像素数量并转换为实际面积：
import numpy as np from rasterio.mask import mask def calculate_water_area(year, boundary): # 加载对应年份数据 img_path = f"JRC_GSW_{year}.tif" with rasterio.open(img_path) as src: # 精确裁剪 clipped, transform = mask(src, boundary.geometry, crop=True) water_mask = (clipped[0] == 1) # 1表示水体 # 计算面积（考虑地球曲率） pixel_area = abs(transform[0] * transform[4]) * 1e6 # 转换为平方米 total_area = np.sum(water_mask) * pixel_area return total_area # 计算1984-2019年序列 years = range(1984, 2020) area_series = [calculate_water_area(y, city_boundary) for y in years]3.3 常见问题解决方案
投影转换问题：当出现"transform doesn't match"错误时，需统一所有数据的CRS：
# 强制重投影示例 with rasterio.open('input.tif') as src: transform, width, height = calculate_default_transform( src.crs, 'EPSG:32650', src.width, src.height, *src.bounds) kwargs = src.meta.copy() kwargs.update({ 'crs': 'EPSG:32650', 'transform': transform, 'width': width, 'height': height }) with rasterio.open('reprojected.tif', 'w', **kwargs) as dst: reproject( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs='EPSG:32650', resampling=Resampling.nearest)4. 动态可视化呈现
4.1 制作变迁动图
使用matplotlib的动画模块生成GIF：
from matplotlib.animation import FuncAnimation import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(10, 8)) def update(frame): ax.clear() year = years[frame] img = rasterio.open(f"JRC_GSW_{year}.tif") show(img, ax=ax, cmap='Blues') ax.set_title(f"{year}年水域分布", fontsize=14) ani = FuncAnimation(fig, update, frames=len(years), interval=500) ani.save('water_change.gif', writer='pillow', dpi=150)4.2 交互式地图展示
结合folium创建可缩放的时间轴地图：
import folium from folium.plugins import TimeSliderChoropleth # 创建底图 m = folium.Map(location=[30.25, 120.2], zoom_start=10) # 准备时间序列数据 style_function = lambda x: { 'fillColor': '#1f77b4' if x['properties']['is_water'] else '#d3d3d3', 'weight': 0.5 } TimeSliderChoropleth( data=geojson_data, styledict=style_dict, overlay=True, control=True ).add_to(m) m.save('interactive_map.html')5. 深度分析技巧
5.1 变化热点检测
通过Change图层识别显著变化区域：
# 计算变化率前10%的区域 change_layer = rasterio.open('JRC_GSW_Change.tif') change_data = change_layer.read(1) threshold = np.percentile(change_data[change_data>0], 90) hotspots = np.where(change_data >= threshold) print(f"发现{len(hotspots[0])}个显著变化热点")5.2 结合OpenStreetMap标注
将分析结果与现实地物对应：
import osmnx as ox # 获取西湖周边建筑数据 west_lake = ox.geometries_from_place('西湖, 杭州', tags={'building':True}) # 叠加水域变化图层 fig, ax = plt.subplots() water_change.plot(ax=ax, alpha=0.5) west_lake.plot(ax=ax, color='red', markersize=5)处理这类地理数据时最常遇到的坑是坐标系统不匹配——全球数据集通常使用WGS84(EPSG:4326)，而本地规划数据可能采用高斯-克吕格投影。建议在流程开始时就统一所有数据的CRS，可以节省大量调试时间。