尧图网站建设 
简介新疆维吾尔自治区水库时空数据集覆盖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_year1965的水库去跟 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:4214BJZ54 地理坐标系再做七参数转换。不过由于历史资料本身精度低很多情况下也可以直接按 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_overrideTrue) # 转换到 CGCS2000 地理坐标 transformer Transformer.from_crs(EPSG:4214, EPSG:4490, always_xyTrue) 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) # 统一到 WGS84EPSG:4326方便与遥感切片叠加 gdf_4326 gdf_raw.to_crs(EPSG:4326)这里要说明allow_overrideTrue只用于.prj缺失、且你明确知道源坐标系的情况。always_xyTrue保证经纬度顺序为 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, encodingutf-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, axis1) # 检查重复编码 dup df_state[df_state.duplicated(reservoir_id, keepFalse)] 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(timendwi.time.dt.month.isin([7, 8])) # 使用 groupby 取年际中值避免云污染造成空值 ndwi_median ndwi_yearly.where(ndwi_yearly 0.5).groupby(time.year).median(dimtime) # 按阈值提取水域多边形 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(yearyear).squeeze().values mask arr 0.2 # 调整阈值 results shapes(arr, maskmask, transformndwi_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, crsEPSG:4326) # 空间关联到具体水库 gdf_reservoirs gpd.read_file(clean/reservoir_entities.gpkg) join gpd.sjoin(gdf_extents, gdf_reservoirs[[reservoir_id, geometry]], howinner, predicateintersects) 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(indexreservoir_id, columnsyear, valuesarea_km2) # 填充缺失值这里用相邻年份均值 pivot pivot.interpolate(axis1, limit6).fillna(0) # 归一化保留振幅差异时用 MinMaxScaler data_scaled TimeSeriesScalerMeanVariance().fit_transform(pivot.values) # 设置聚类数k4可先通过 silhouette score 确定 model TimeSeriesKMeans(n_clusters4, metricdtw, max_iter50, random_state0) labels model.fit_predict(data_scaled) # 输出每个水库的聚类标签和中心趋势 cluster_df pd.DataFrame({reservoir_id: pivot.index, cluster: labels}) # 打印每类的数量 print(cluster_df.groupby(cluster).size())metricdtw不适合长短不一的序列但这里已经是年对齐的等长序列。更常用的做法是先用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放在同一个事务里避免统计信息滞后导致执行计划走偏。本文还有配套的精品资源点击获取 构建水库时空数据集：四层模型与时空SQL分析实战
简介新疆维吾尔自治区水库时空数据集覆盖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_year1965的水库去跟 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:4214BJZ54 地理坐标系再做七参数转换。不过由于历史资料本身精度低很多情况下也可以直接按 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_overrideTrue) # 转换到 CGCS2000 地理坐标 transformer Transformer.from_crs(EPSG:4214, EPSG:4490, always_xyTrue) 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) # 统一到 WGS84EPSG:4326方便与遥感切片叠加 gdf_4326 gdf_raw.to_crs(EPSG:4326)这里要说明allow_overrideTrue只用于.prj缺失、且你明确知道源坐标系的情况。always_xyTrue保证经纬度顺序为 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, encodingutf-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, axis1) # 检查重复编码 dup df_state[df_state.duplicated(reservoir_id, keepFalse)] 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(timendwi.time.dt.month.isin([7, 8])) # 使用 groupby 取年际中值避免云污染造成空值 ndwi_median ndwi_yearly.where(ndwi_yearly 0.5).groupby(time.year).median(dimtime) # 按阈值提取水域多边形 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(yearyear).squeeze().values mask arr 0.2 # 调整阈值 results shapes(arr, maskmask, transformndwi_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, crsEPSG:4326) # 空间关联到具体水库 gdf_reservoirs gpd.read_file(clean/reservoir_entities.gpkg) join gpd.sjoin(gdf_extents, gdf_reservoirs[[reservoir_id, geometry]], howinner, predicateintersects) 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(indexreservoir_id, columnsyear, valuesarea_km2) # 填充缺失值这里用相邻年份均值 pivot pivot.interpolate(axis1, limit6).fillna(0) # 归一化保留振幅差异时用 MinMaxScaler data_scaled TimeSeriesScalerMeanVariance().fit_transform(pivot.values) # 设置聚类数k4可先通过 silhouette score 确定 model TimeSeriesKMeans(n_clusters4, metricdtw, max_iter50, random_state0) labels model.fit_predict(data_scaled) # 输出每个水库的聚类标签和中心趋势 cluster_df pd.DataFrame({reservoir_id: pivot.index, cluster: labels}) # 打印每类的数量 print(cluster_df.groupby(cluster).size())metricdtw不适合长短不一的序列但这里已经是年对齐的等长序列。更常用的做法是先用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放在同一个事务里避免统计信息滞后导致执行计划走偏。本文还有配套的精品资源点击获取