用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的水体变化程度Extent35年间曾出现过水体的所有位置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 python3.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(edgecolorred, facecolornone) ctx.add_basemap(ax, crscity_boundary.crs, sourcectx.providers.Stamen.Terrain)3.2 年度水域面积计算通过遍历每年数据统计水域像素数量并转换为实际面积import numpy as np from rasterio.mask import mask def calculate_water_area(year, boundary): # 加载对应年份数据 img_path fJRC_GSW_{year}.tif with rasterio.open(img_path) as src: # 精确裁剪 clipped, transform mask(src, boundary.geometry, cropTrue) 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 doesnt 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( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsEPSG:32650, resamplingResampling.nearest)4. 动态可视化呈现4.1 制作变迁动图使用matplotlib的动画模块生成GIFfrom 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(fJRC_GSW_{year}.tif) show(img, axax, cmapBlues) ax.set_title(f{year}年水域分布, fontsize14) ani FuncAnimation(fig, update, frameslen(years), interval500) ani.save(water_change.gif, writerpillow, dpi150)4.2 交互式地图展示结合folium创建可缩放的时间轴地图import folium from folium.plugins import TimeSliderChoropleth # 创建底图 m folium.Map(location[30.25, 120.2], zoom_start10) # 准备时间序列数据 style_function lambda x: { fillColor: #1f77b4 if x[properties][is_water] else #d3d3d3, weight: 0.5 } TimeSliderChoropleth( datageojson_data, styledictstyle_dict, overlayTrue, controlTrue ).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_data0], 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(axax, alpha0.5) west_lake.plot(axax, colorred, markersize5)处理这类地理数据时最常遇到的坑是坐标系统不匹配——全球数据集通常使用WGS84(EPSG:4326)而本地规划数据可能采用高斯-克吕格投影。建议在流程开始时就统一所有数据的CRS可以节省大量调试时间。