用JRC全球地表水数据，5分钟搞定你所在城市的水体变迁分析（附Python代码）

用JRC全球地表水数据5分钟分析城市水体变迁（附Python实战代码）

每次回老家总听长辈念叨："小时候村口那条河比现在宽多了"、"东边的水库是后来才修的"。作为地理信息爱好者，我总想用数据验证这些记忆。JRC全球地表水数据集就像一台时空望远镜，让我们能直观看到1984年以来每片水域的变迁轨迹。今天分享的这套方法，用不到10行Python代码就能生成专业级水体变化分析图。

1. 准备工作：数据获取与环境配置

1.1 数据集快速下载指南

JRC全球地表水数据集包含7个关键图层，我们重点关注：

• extent（水体存在标识）：二进制标记1984-2019年间是否出现过水体
• change（变化强度）：用0-254的数值量化水体增减程度

通过Google Earth Engine直接获取数据最便捷，若需本地处理可下载GeoTIFF格式分幅数据。推荐使用以下命令快速获取目标区域数据：

import ee ee.Initialize() # 定义杭州区域边界 hangzhou = ee.Geometry.Rectangle([119.2, 29.4, 120.8, 30.5]) # 获取2019年水体分布数据 water_2019 = ee.ImageCollection("JRC/GSW1_2/MonthlyHistory")\ .filterBounds(hangzhou)\ .mosaic()\ .clip(hangzhou)

1.2 Python环境速配方案

新建conda环境并安装核心库：

conda create -n water_analysis python=3.8 conda activate water_analysis pip install geopandas rasterio matplotlib earthpy

注意：若使用本地下载的数据文件，需确保GDAL库版本≥3.0，可通过conda install -c conda-forge gdal安装

2. 数据解析：理解JRC图层编码规则

2.1 extent图层实战解读

这个二值图层就像水体"存在证明"，每个30m×30m的网格只有两种状态：

用rasterio快速查看数据分布：

import rasterio from rasterio.plot import show with rasterio.open('water_extent.tif') as src: print(f"水体覆盖率：{src.read(1).mean()*100:.2f}%") show(src, title='历史水体分布')

2.2 change图层深度解码

变化强度图层采用特殊编码体系，需要转换才能得到直观百分比：

import numpy as np def decode_change(arr): """将原始编码转换为变化百分比""" result = np.zeros_like(arr, dtype=np.float32) result[arr==0] = -1 # 100%减少 result[arr==100] = 0 # 无变化 result[arr==200] = 1 # 100%增加 return result change_array = rasterio.open('change.tif').read(1) percent_change = decode_change(change_array)

3. 城市水体变迁分析实战

3.1 变化热点区域定位

结合geopandas进行空间统计，快速找出变化最显著的区域：

import geopandas as gpd from shapely.geometry import shape # 生成变化强度等值线 contours = gpd.GeoDataFrame.from_features( earthpy.spatial.contour_array( percent_change, interval=0.2, transform=src.transform ) ) # 筛选显著变化区域 hotspots = contours[contours['value'].abs() > 0.5] print(f"发现{len(hotspots)}个显著变化区域")

3.2 时空变化趋势可视化

使用matplotlib制作专业级分析图：

import matplotlib.pyplot as plt fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6)) # 左侧绘制历史水体分布 extent = src.read(1) ax1.imshow(extent, cmap='Blues') ax1.set_title('1984-2019水体存在记录') # 右侧绘制变化强度 im = ax2.imshow(percent_change, cmap='coolwarm', vmin=-1, vmax=1) plt.colorbar(im, ax=ax2, label='变化强度') ax2.set_title('水体变化趋势') plt.savefig('water_change_analysis.png', dpi=300)

4. 高级技巧：多时段对比分析

4.1 年代际变化统计

将1984-1999与2000-2019两个时段数据对比：

# 计算不同年代水体面积 def calc_decade_area(year_range): decade_extent = (extent & (year_mask(year_range))).astype(int) return decade_extent.sum() * 900 / 1e6 # 转换为平方公里 areas = { "1984-1999": calc_decade_area((1984,1999)), "2000-2019": calc_decade_area((2000,2019)) } print(pd.Series(areas).to_markdown())

输出示例：

4.2 季节性水体识别

结合seasonality图层分析水体稳定性：

seasonal = rasterio.open('seasonality.tif').read(1) permanent_water = (seasonal == 12) # 全年存在水体 seasonal_water = (seasonal > 0) & (seasonal < 12) plt.figure(figsize=(10,6)) plt.imshow(permanent_water, cmap='Blues', alpha=0.5) plt.imshow(seasonal_water, cmap='Oranges', alpha=0.5) plt.legend(['永久水体','季节性水体'])

5. 常见问题解决方案

5.1 坐标系转换问题

当数据坐标参考系与行政边界不匹配时：

from rasterio.warp import calculate_default_transform, reproject def reproject_raster(src_file, dst_crs='EPSG:4326'): """重投影栅格数据""" with rasterio.open(src_file) as src: transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs = src.meta.copy() kwargs.update({ 'crs': dst_crs, '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=dst_crs)

5.2 大数据量处理技巧

对于市级以上区域，建议使用分块处理：

from rasterio.windows import Window def process_by_chunk(filename, chunk_size=1024): """分块处理大文件""" with rasterio.open(filename) as src: for ji, window in src.block_windows(): chunk = src.read(window=window) # 在此处添加处理逻辑 yield window, processed_chunk

记得第一次跑通整个流程时，我对着生成的水体变化图发了半天呆——原来小区旁边那个公园人工湖在2008年前根本不存在。这种用数据验证直觉的过程，正是地理分析的魅力所在。建议初学者先从家乡熟悉的区域开始分析，当数据与记忆产生碰撞时，往往会有意外发现。