简介：新疆维吾尔自治区水库时空数据集覆盖1942—2022年，面向地理信息、水利规划与区域发展研究者，可支撑长时序水库演变分析、流域对比及规划选址决策。基于2022年Sentinel-2影像及历史水库记录，提取全区面积大于0.001平方千米的水库空间分布，并整理水库名称、经纬度、平均海拔、建成年份、总库容及最大面积、流域等属性字段。资源包共9个文件，以.shp矢量格式及配套文件（.shx、.dbf、.prj等）为主，另有.xlsx属性表，压缩后大小5.87MB，可直接在ArcGIS或QGIS中加载，便于开展空间查询与制图，已有178人学习下载。数据按水电标准划分大、中、小型水库，截至2022年共收录804座，总库容24.16立方千米，其中大型37座、中型175座、小型592座，平原水库461座、山地水库343座。同时提供1942—2022年逐期空间分布成果，可直观反映新疆水库建设时空动态，为区域水资源管理与科学研究提供可靠数据基础。
1. 为什么要把 80 年的水库绑进一个时空数据集
做水利信息化或者灾备规划的人大概都有同感：新疆的水库台账散得厉害，1942 年建成的老库可能只留在纸质竣工图里，2000 年后的新库才有完整的位置、库容和调度记录。想回答“一个流域 80 年里到底新增了多少座库容大于 1000 万方的水库”“哪些老库的淹没范围变了”这类问题，光靠一张现状矢量图远远不够——你需要的是把时间维度显式地建模进去。水库时空数据集，本质上是把每一次“水库出现、扩容、废弃、位置修正”的事件，连同它的几何边界和属性字段，存储成可查询、可计算、可可视化的时序地理实体。本文从数据结构、ETL 流程、时空分析三个方面展开，适合刚接触时空数据建模的 GIS 工程师，也适合需要处理长序列水利数据的后端开发。
这套数据集的价值不在“多一张图”，而在于它把离散的历史资料统一到一个时间坐标系里，让你能用时空 SQL 和 Python 直接做趋势统计、变化检测和蓄水能力回溯。下面从数据怎么组织开始讲。
2. 水库时空数据集的四层模型与目录设计
2.1 时空数据集的“四层”分法
我处理过的水利时空数据，通常不会只存一张大宽表，而是拆成四层，每一层解决一个类型的问题。第一层是基础地理层，包括流域边界、河流中心线、行政区划，这些在时间上基本不变，用普通矢量数据存即可。第二层是水库实体层，每个水库有一个全局唯一标识（水库编码），描述它的名称、所在河流、建成时间、设计库容、坝型、功能属性，这部分是“实体主档”。第三层是时间状态层，记录每个水库在某个时间点或时间段的状态变更，比如 1965 年扩建后库容从 500 万方变成 1200 万方，或者 1988 年因为淤积导致有效库容减少。第四层是空间形态层，存水库的边界多边形或坝址点，但空间形状本身也可能随时间变化——老水库在枯水期和丰水期的水域面积差异很大，所以空间层还要区分“坝址点位”和“水域范围”两种几何。
这个分法的好处是：查询“某年某流域有几个大型水库”时，不需要去所有时间切片里做空间叠加，只要在实体层过滤状态层；而要做“水域面积变化”时，直接查空间形态层的时间序列。下面给出一个我常用的目录结构。
reservoir_spatiotemporal/ ├── raw/ # 原始资料，不可修改 │ ├── historical_docs/ # 纸质档案扫描件、PDF │ ├── survey_shp/ # 不同年代测绘的shapefile │ └── satellite/ # Landsat等影像，按条带组织 ├── clean/ # 清洗后的中间数据 │ ├── reservoir_entities.gpkg │ ├── reservoir_state.csv │ └── water_extents/ # 按提取年份分层的多边形 ├── feature/ # 分析特征 │ ├── reservoir_timeline.parquet │ └── metrics/ # 库容、面积、蓄水量等指标 └── export/ # 对外发布的GeoJSON或SQLiteraw 目录只读，clean 目录记录清洗脚本的版本，feature 目录是给下游分析和机器学习用的宽表。这样当别人问“你这个面积数据是哪一年的”，你能快速回溯到卫星影像和提取脚本。
2.2 时间属性的表示：三种模式要分清
水库时空数据集里最容易出错的是时间表示。我看到很多项目把“建成时间”当成唯一时间字段，然后所有分析都基于这个字段，这就丢了“空间观测时间”。我一般把时间分为事件时间、观测时间和状态有效期三类。事件时间是水库发生某次变更（如竣工、除险加固）的时点，用valid_from和valid_to两个字段表示状态有效期；观测时间是某次遥感影像或实地测绘实际采样的日期，用survey_date字段表示；分析参考时间是你要切片查询的时刻，这通常不是存储字段，而是你在查询条件里传的参数。
这三种时间混在一起会导致查询结果失真。举个例子：如果你用build_year=1965的水库去跟 2020 年的影像做叠加，得到的水域面积其实是 2020 年的现状，而不是 1965 年该水库建成时的面积。正确的做法是单独建一个“水库状态历史表”，每一行是一个版本记录。
-- PostgreSQL + PostGIS 示例 CREATE TABLE reservoir_state_history ( reservoir_id varchar(32) NOT NULL, valid_from date NOT NULL, -- 状态生效日期 valid_to date, -- 状态失效日期，NULL表示当前 capacity_10000m3 numeric(12,2), storage_10000m3 numeric(12,2), dead_capacity_10000m3 numeric(12,2), status varchar(16), -- 运行/在建/废弃/改建 PRIMARY KEY (reservoir_id, valid_from) ); CREATE TABLE reservoir_extent_timeline ( extent_id serial PRIMARY KEY, reservoir_id varchar(32) REFERENCES reservoir_state_history(reservoir_id), survey_date date NOT NULL, -- 实际观测日期 water_area_km2 numeric(10,4), geom geometry(Polygon, 4326) -- 水域边界 );注意valid_to为 NULL 代表当前有效，查询时用(valid_from <= :query_date AND (valid_to IS NULL OR valid_to > :query_date))。如果你用valid_to存一个大日期如9999-12-31，也可以，但要注意部分统计函数会把非 NULL 值算进去，容易出错。
2.3 数据源的时空分辨率统一
新疆的水库数据来源跨度极大，1942 年的资料可能是 1:5 万地形图上标出的一个点，1990 年代有 TM 影像解译的边界，2020 年后有高分二号影像和无人机测绘。把这些数据放一起，必须明确各自的时空分辨率，否则叠加后会出现大量拓扑错误。我的做法是给每一条空间数据加两个元数据字段：source_scale和positional_accuracy_m。前者记录比例尺或 GSD，后者记录估计的位置精度（单位米）。历史点位数据我一般给 100–500 米，现代影像解译给 10–30 米。后续做空间分析时，如果两个数据源的位置精度相差超过阈值，我会做“降精度对齐”——把高精度几何缓冲到低精度半径再参与计算，避免因为 1 米的边界让一个 1942 年的点和 2022 年的线“擦边不交”。
3. 从历史资料到时空入库的 ETL 实战
3.1 历史矢量资料的地理配准和编码修复
新疆地区 1950 年代前后的测绘成果多采用北京 1954 坐标系（也就是常说的 BJZ54），而现代影像和 GNSS 数据用的是 CGCS2000 或 WGS84。处理流程第一步是统一坐标系。我推荐用 PROJ 库操作，而不是在 ArcGIS 里手动转一遍。因为矢量数据量大，常见做法是先用 GDAL 检查原始坐标系描述是否完整。
# 检查 shapefile 的 .prj 文件是否含完整参数 ogrinfo -al -so raw/survey_shp/xinjiang_reservoir_1958.shp | grep -i "CRS" # 如果没有投影信息，使用 gdal.SetGeometry 时需要指定源坐标系Python 脚本中，用 GeoPandas 读取后统一转成 EPSG:4326。注意，如果原始数据是老的克拉索夫斯基椭球下的 BJZ54，要使用明确转换方法，不要直接用to_crs("EPSG:4326")，因为可能默认套用 WGS84 转换导致水平误差达到几十米。正确做法是先定义为EPSG:4214（BJZ54 地理坐标系），再做七参数转换。不过由于历史资料本身精度低，很多情况下也可以直接按 CF 网格做近似转换，但要记录在元数据里。
import geopandas as gpd from pyproj import CRS, Transformer # 读取原始数据，强制指定源坐标系为 BJZ54 地理坐标 gdf_raw = gpd.read_file("raw/survey_shp/xinjiang_reservoir_1958.shp") # 如果 .prj 缺失，需要手动设置： gdf_raw = gdf_raw.set_crs("EPSG:4214", allow_override=True) # 转换到 CGCS2000 地理坐标 transformer = Transformer.from_crs("EPSG:4214", "EPSG:4490", always_xy=True) gdf_raw["geometry"] = gdf_raw.geometry.map( lambda geom: transform_geom(geom, transformer) ) def transform_geom(geom, transformer): from shapely.ops import transform import shapely if geom.is_empty: return geom return transform(lambda x, y: transformer.transform(x, y), geom) # 统一到 WGS84（EPSG:4326）方便与遥感切片叠加 gdf_4326 = gdf_raw.to_crs("EPSG:4326")这里要说明：allow_override=True只用于.prj缺失、且你明确知道源坐标系的情况。always_xy=True保证经纬度顺序为 x=经度、y=纬度。之后用to_crs("EPSG:4326")把 CGCS2000 转 WGS84，两者在新疆地区差异不到 1 米，对历史数据来说可忽略。
3.2 属性字段的清洗与水库编码生成
原始台账里的地名经常不一致，比如“肯斯瓦特水库”可能在不同年代写成“克斯瓦特水库”。解决方法是建立别名映射表，并用编码中心发号的模式生成 12 位水库编码。编码结构为：6 位流域编码 + 2 位建成年份后两位 + 4 位序号。这样从编码里就能直接看出流域和大致年代。
import pandas as pd import hashlib # 读取多期合并的属性表 df_state = pd.read_csv("clean/tmp_reservoir_all.csv") # 使用别名映射表统一名称 alias_map = pd.read_csv("ref/reservoir_alias.csv", encoding="utf-8") name_to_canonical = dict(zip(alias_map["alias"], alias_map["canonical_name"])) df_state["canonical_name"] = df_state["reservoir_name"].map( lambda x: name_to_canonical.get(str(x).strip(), str(x).strip()) ) # 生成水库编码：流域代码 + 建成年份 + 四角坐标哈希尾数 def make_reservoir_id(row): basin_code = row["basin_code"] # 如 'II' 代表额尔齐斯河流域 yr = str(row["built_year"])[-2:] if pd.notna(row["built_year"]) else "00" # 用坐标字符串生成稳定哈希作为序号 coord_key = f"{row['lon']:.4f},{row['lat']:.4f}" seq = hashlib.md5(coord_key.encode()).hexdigest()[:4].upper() return f"{basin_code}{yr}{seq}" df_state["reservoir_id"] = df_state.apply(make_reservoir_id, axis=1) # 检查重复编码 dup = df_state[df_state.duplicated("reservoir_id", keep=False)] if len(dup): print("重复编码记录：", len(dup)) # 处理逻辑：追加序号或人工审核这里有个容易忽略的坑：同一座水库如果建成年份推算有误，生成的编码就会变化，导致后续历史状态无法关联。我的习惯是先生成一个“候选实体表”，用名称和坐标做空间匹配（距离小于 500 米就认为是同一实体），再人工确认。
3.3 栅格时间序列的整理与入库
水库空间形态的变化通常从卫星影像提取水域面积，这里以 Landsat 系列为例。用xarray把不同期影像组织成带时间的数组，然后基于 NDWI 提取水体。注意，新疆的山区水库在 4–5 月融雪期和 8–9 月汛期水面差异很大，如果数据集时间跨度大，却只选了一期的影像，那这个面积不能代表该年份的典型状态。我建议在提取脚本里加一个时间筛选参数。
import xarray as xr import rioxarray from datetime import datetime # 构建一个包含 time 维的 zarr 数据集（示意路径） da_nir = xr.open_zarr("clean/satellite/nir.zarr")["band1"] da_g = xr.open_zarr("clean/satellite/green.zarr")["band1"] # 计算 NDWI，并提取每年 7-8 月的中位数影像作为“丰水期代表” ndwi = (da_g - da_nir) / (da_g + da_nir) ndwi_yearly = ndwi.sel(time=ndwi.time.dt.month.isin([7, 8])) # 使用 groupby 取年际中值，避免云污染造成空值 ndwi_median = ndwi_yearly.where(ndwi_yearly < 0.5).groupby("time.year").median(dim="time") # 按阈值提取水域多边形 from rasterio.features import shapes import geopandas as gpd from shapely.geometry import shape # 对每一年的二维切片进行矢量化（简化示意） all_extents = [] for year in ndwi_median.year.values: arr = ndwi_median.sel(year=year).squeeze().values mask = arr > 0.2 # 调整阈值 results = shapes(arr, mask=mask, transform=ndwi_median.rio.transform()) geoms = [shape(g) for g, v in results if v == 1] # 过滤面积大于 0.1 km² 的水体 for geom in geoms: if geom.area * 111 * 111 > 0.1: # 粗略换算 all_extents.append({ "year": int(year), "area_km2": geom.area * 111 * 111, "geometry": geom, }) gdf_extents = gpd.GeoDataFrame(all_extents, crs="EPSG:4326") # 空间关联到具体水库 gdf_reservoirs = gpd.read_file("clean/reservoir_entities.gpkg") join = gpd.sjoin(gdf_extents, gdf_reservoirs[["reservoir_id", "geometry"]], how="inner", predicate="intersects") join = join[join.area_km2 > 0.1]ndwi_yearly.where(ndwi_yearly < 0.5)是为了去除裸地和盐碱地的高值噪声。这个阈值和mask = arr > 0.2需要根据具体影像调整，我在新疆南部用 Landsat-8 时通常把阈值固定在 0.15，但经过地形阴影校正后可能要改为 0.25。这里给的是最直接的做法，实际项目中建议先对几个典型水库做目视比较。
4. 用时空 SQL 与 Python 做水库变化分析
4.1 按“有效日期”查询水库状态历史
这是时空数据集最核心的查询模式。存在reservoir_state_history表后，任何时间点的“水库快照”都可以用一个标准 SQL 取出来。
-- 查询1980年12月31日时，各流域的大型水库数量（库容≥1亿m³） WITH snap AS ( SELECT DISTINCT ON (reservoir_id) reservoir_id, capacity_10000m3, valid_from, valid_to FROM reservoir_state_history WHERE valid_from <= DATE '1980-12-31' AND (valid_to IS NULL OR valid_to > DATE '1980-12-31') ORDER BY reservoir_id, valid_from DESC, valid_to NULLS LAST ) SELECT substr(reservoir_id, 1, 6) AS basin_code, count(*) AS reservoir_count, sum(capacity_10000m3) AS total_capacity FROM snap WHERE capacity_10000m3 >= 10000 -- 库容单位：万m³，1亿m³即10000万m³ GROUP BY substr(reservoir_id, 1, 6) ORDER BY total_capacity DESC;这里的关键是DISTINCT ON (reservoir_id)，它保证每个水库只取一条符合条件的记录。注意valid_to IS NULL表示记录至今有效。substr(reservoir_id, 1, 6)用前缀提取流域编码。如果你用的是日期类型，不要用字符串比较，否则索引会失效。
4.2 水域面积时序的自相关修正
从reservoir_extent_timeline做面积变化分析时，不能直接算年际差值，因为不同年份的观测日期不一样，丰枯水期造成的面积波动可能远大于真实扩张。一个常见做法是计算“同月年均值”，即把每年 7–8 月的观测取平均作为“高水位面积”，再用 10–11 月作为“低水位面积”。在 SQL 里可以这样：
SELECT reservoir_id, date_part('year', survey_date) AS year, count(*) AS obs_count, avg(water_area_km2) AS avg_high_area FROM reservoir_extent_timeline WHERE date_part('month', survey_date) BETWEEN 7 AND 8 GROUP BY reservoir_id, date_part('year', survey_date) ORDER BY reservoir_id, year;如果某一年多云导致缺测，我一般不做线性插值，而是用前后两年的同月均值填充，并标记is_interpolated字段。否则后续做趋势检验时，插值造成的假平滑会让你误判“显著变化”。
4.3 时序聚类找出异常演变模式
水库时空数据集中最有价值的分析之一是把所有水库的库容/面积时间序列做聚类，识别出“持续扩容型”“淤积萎缩型”“周期性波动型”等模式。下面给一个用 tslearn 做形状聚类的脚本框架。
import pandas as pd import numpy as np from tslearn.clustering import TimeSeriesKMeans from tslearn.preprocessing import TimeSeriesScalerMeanVariance # 读取特征宽表：索引是水库编码，列是年份或观测期 df_wide = pd.read_parquet("feature/reservoir_area_long.parquet") # 按水库ID和年份整理成矩阵 pivot = df_wide.pivot(index="reservoir_id", columns="year", values="area_km2") # 填充缺失值，这里用相邻年份均值 pivot = pivot.interpolate(axis=1, limit=6).fillna(0) # 归一化，保留振幅差异时用 MinMaxScaler data_scaled = TimeSeriesScalerMeanVariance().fit_transform(pivot.values) # 设置聚类数k=4，可先通过 silhouette score 确定 model = TimeSeriesKMeans(n_clusters=4, metric="dtw", max_iter=50, random_state=0) labels = model.fit_predict(data_scaled) # 输出每个水库的聚类标签和中心趋势 cluster_df = pd.DataFrame({"reservoir_id": pivot.index, "cluster": labels}) # 打印每类的数量 print(cluster_df.groupby("cluster").size())metric="dtw"不适合长短不一的序列，但这里已经是年对齐的等长序列。更常用的做法是先用min-max归一化再做欧氏距离聚类，因为水库面积的绝对数值对工程判断很重要，比如一个 5 平方公里的水库面积波动 1 平方公里和一个 0.5 平方公里的水库波动 0.2 平方公里不能直接比，“时序缩放”会抹掉这种量级差异，所以这里我用了TimeSeriesScalerMeanVariance，保留均值信息。
4.4 空间相邻水库的联合变化检测
如果想知道两个相邻水库是否出现“你涨我消”的博弈，可以用空间连接加时间相关性计算。PostGIS 里可以先做ST_DateLine空间缓冲区连接，再在 Python 里按对计算 Pearson 相关系数。
# 计算每个水库与其他水库的最小距离（单位：度） psql -d water_db -c "SELECT a.reservoir_id, b.reservoir_id, ST_Distance(a.geom, b.geom) AS dist_deg FROM reservoir_extent_timeline a JOIN reservoir_extent_timeline b ON a.reservoir_id < b.reservoir_id WHERE a.obs_date = '2020-08-01' AND b.obs_date = '2020-08-01' AND ST_DWithin(a.geom::geography, b.geom::geography, 30000) ORDER BY dist_deg LIMIT 100;"ST_DWithin用::geography转成地理坐标后再按米计算，比直接算平面度准确，但在新疆这种高纬度地区，30 公里的平面度误差很小，所以也可以用ST_Distance先做粗筛。
5. 数据质量验证与长序列查询的索引优化技巧
5.1 三条常用质量校验规则
首先检查时间闭合性。状态历史表中每个水库的valid_from和valid_to应该首尾相接，不能出现断裂或重叠。可以写一个 SQL 检查重叠：
SELECT a.reservoir_id, a.valid_from, a.valid_to, b.valid_from, b.valid_to FROM reservoir_state_history a JOIN reservoir_state_history b ON a.reservoir_id = b.reservoir_id AND a.valid_from < b.valid_from AND a.valid_to > b.valid_from AND a.valid_to < b.valid_to;这条JOIN如果返回任何行，说明存在交叉覆盖。修复方式是把后一条的valid_from调整为前一条的valid_to。
其次检查空间拓扑。提取的水域多边形不能互相重叠，尤其是同一水库的相邻年份多边形，正常情况应该大致包含或相交，但不能出现一个多边形完全吞掉另一个后又变小的诡异情况。用 PostGISST_Intersects配合ST_Area计算重叠面积占比。
最后检查属性逻辑。库容和面积应大致满足幂函数关系，如果某个水库面积很大但库容很小，很可能是面积提取包含了浅滩或盐沼。做回归时剔除异常值。
5.2 时空索引的构建顺序
长序列查询最怕全表扫描，PostGIS 中GiST索引需要同时覆盖空间和时间字段。我常用的建索引 SQL 是：
CREATE INDEX idx_extent_space ON reservoir_extent_timeline USING gist (geom); CREATE INDEX idx_extent_time ON reservoir_extent_timeline USING btree (survey_date); CREATE INDEX idx_extent_comp ON reservoir_extent_timeline USING gist (geom, survey_date);第一条加速空间过滤，第二条加速时间范围过滤，第三条是组合索引，用于类似“某空间范围内某时间段的水库变化”查询。注意GiST组合索引的查询顺序要尽量与索引字段顺序一致，即查询条件先匹配geom，再匹配survey_date，所以写 SQL 时WHERE geom && :bbox AND survey_date BETWEEN ...是比较理想的。
5.3 用物化视图缓存年度快照
如果查询频率很高的“逐十年快照统计”已经固定，我建议建一个物化视图定期刷新，避免每次都对整个历史表做DISTINCT ON。例如：
CREATE MATERIALIZED VIEW reservoir_snapshot_1990 AS SELECT DISTINCT ON (reservoir_id) reservoir_id, capacity_10000m3, status FROM reservoir_state_history WHERE valid_from <= DATE '1990-12-31' AND (valid_to IS NULL OR valid_to > DATE '1990-12-31') ORDER BY reservoir_id, valid_from DESC;这个物化视图可以按十年打几个固定快照，也可以只建当前年份的。刷新时机选在数据完成增量 ETL 之后。这样你日常做探索分析时，就不需要反复扫描历史表了。一个我习惯的小技巧：把物化视图的刷新和ANALYZE放在同一个事务里，避免统计信息滞后导致执行计划走偏。
本文还有配套的精品资源，点击获取
