ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

1989-2020年中国滨海滩涂湿地遥感数据集处理与面积统计指南

1989-2020年中国滨海滩涂湿地遥感数据集处理与面积统计指南 简介针对中国北纬18度以北沿海滩涂湿地这份数据集收录了1989—2020年间32个年份的分布时空变化信息适合生态学、地理学、环境科学及城市规划领域的研究人员用于分析湿地面积增减、岸线迁移与退化恢复趋势。压缩包内共256个文件为完整shapefile格式由shp、dbf、prj、sbn、sbx、shx、cpg、xml等32组配套文件构成可记录各年份湿地的位置、面积、类型等属性支持ArcGIS、QGIS直接读取包体大小约201MB。目前已有342人学习浏览。通过逐年对比可识别人类活动与气候变化影响下的湿地热点区域为湿地保护、海岸带管理和政策制定提供科学依据是理解近30年滨海湿地演变的重要基础数据。1. 18°N以北的33年滩涂影像先弄懂边界和分类再谈用数据“18°N以北中国滨海滩涂湿地分布数据集1989-2020.rar”这个标题拆开看信息密度很高。“18°N”从空间上划定了数据覆盖范围的南界从中国海岸带地理来看辽宁鸭绿江口到广西北仑河口的整条大陆海岸带都在范围内海南岛南部一小段和南海诸岛被排除在外“1989-2020”跨越了Landsat 5后期、Landsat 7中期和Landsat 8后期的完整遥感时序“滨海滩涂湿地”是研究对象也就是潮间带反复出露和淹没的滩面及其上覆植被。把这三件事想清楚再解压后续才能避免面积统计上的方向性错误。这类数据集最常见的用途是海岸带围填海历史溯源、湿地保护修复后的效果评估以及滩涂资源年度储量核算。使用者分两类一类有遥感基础想直接拿分类结果做二次分析另一类没有遥感背景只需要滩涂分布图做空间叠加。两类人对同一份数据的风险点不同前者容易忽略潮位差异造成的误差边界后者容易把Landsat 7条带缺失和云遮挡产生的空洞当成真实变化。所以下文从数据构成、读取方式、异常排查到面积计算完整走一遍。动手之前有一个建议拿到压缩包先别急着解压先找包内的元数据说明。多数版本会附带txt或word文档记录坐标系、分类体系、使用影像的传感器和潮位筛选条件。这个文件决定后续所有处理参数即使只做最简单的面积统计也要保留原文备查。提示先读元数据再决定用什么工具打开。跳过这一步后面所有面积统计都可能建立在错误的坐标系假设上。2. 滩涂湿地分类口径与1989-2020时序数据的构建逻辑2.1 数据里的滩涂是什么潮间带出露滩面与植被覆盖类型滨海滩涂湿地在遥感分类语境下通常指潮间带上周期性出露的裸滩以及生长在裸滩上的盐沼植被和红树林空间范围大体位于平均高潮线与平均低潮线之间。不同版本的数据集对“湿地”的外沿定义差异很大有的会把潮上带的盐碱地和人工养殖塘划进湿地范围有的只保留自然潮间带滩地。如果你是要做年际面积对比一旦分类口径在时间序列中间发生调整面积对比就失去了比较基础。所以开箱第一步是查看属性字段里TYPE、CLASS_ID这类分类列把所有类别取值打印出来与元数据图例逐一核对。常见分类可以归纳成四类后面所有统计都建议沿用这四类口径裸滩以粉砂、淤泥为主的出露滩面无植被覆盖。盐沼/草滩覆盖碱蓬、芦苇、互花米草等植被的潮间带滩面。红树林分布在华南沿海潮间带的木本植物群落单独归为一类。养殖塘/水田部分版本把人工湿地划进来统计时必须与自然滩涂分开。这四个类别的编码在数据集中通常以整数存储而且不同年份之间可能存在编码不一致。处理年际数据时第一步就是做编码对齐把同一类别在不同年份里的不同编码映射到统一值否则后面分组统计会出现类别数量翻倍的可笑结果。2.2 为什么是1989到202030米分辨率的Landsat时间跨度从1989年到2020年连续32年公开可用的中等分辨率遥感数据里只有Landsat系列能覆盖整段跨度。Landsat从1984年至今持续提供30米分辨率的对地观测重访周期16天南海和东海沿海区域每年都有十几到二十几次有效过境。30米分辨率对滩涂这种大面积连片的地物够用但对宽度只有十几米的小型潮沟会存在明显的混合像元问题边界提取结果会倾向于变粗。2000年之前的高分辨率影像要么不连续、要么获取成本高所以覆盖到这么长时间序列的滩涂分布产品数据源几乎都以Landsat为主。时间段主要可用传感器分辨率对滩涂提取影响最大的问题1989-1999Landsat 5 TM30 m云量覆盖导致年度可用影像少2000-2012Landsat 7 ETM30 m2003年5月后SLC故障条带数据空洞2013-2020Landsat 8 OLI30 m辐射分辨率提升水体边界提取更稳定2003年到2012年是长时序遥感最棘手的时段Landsat 7的扫描线校正器失效后影像边缘出现楔形条带每个条带宽度约12个像元。生产这套数据时通常会采用多期影像合成或邻近年份填充来处理所以如果你在某一期成果上看到条带状的锯齿边界不要直接认定是自然变化应到元数据里找该年份的质量说明字段确认是否做过SLC填补。滩涂提取还有一个容易被忽略的关键环节潮位筛选。Landsat过境是固定时刻传感器只记录当时滩面是否出露。生产方要做的是从当年所有有效影像中挑出低潮位、少云量的一景或几景再结合目视检查提取滩涂外边界。潮位差异在沿海强潮河口可能带来2-3公里的滩涂宽度误差因此年际比较时如果某一年数据比相邻年份明显偏小先怀疑当年筛选时的潮位偏高而不是立刻判定为实际面积减小。2.3 RAR包内部常见的目录组织和文件结构RAR压缩包的体积通常在几百MB到2GB之间内部按年份建目录是最常见的组织方式。这种结构方便逐年抽取也能避免一次性加载全部数据时电脑内存被撑爆。解压后的典型目录长这样1989_2020_coastal_wetland/ ├── 1989/ │ ├── wetland_1989.shp │ ├── wetland_1989.shx │ ├── wetland_1989.dbf │ └── quality_1989.txt ├── 1990/ │ └── ... ├── ... ├── 2020/ ├── wetland_type_dict.xlsx └── metadata.pdf如果提供的是栅格成果目录会改为每年度一个tif文件配合颜色映射表或QGIS样式文件。无论哪种形式文件夹名称里的年份必须和文件名里标注的年份一致不一致的情况虽然少一旦遇到会直接影响后面批量读取时的年份提取逻辑值得单独检查。2.4 坐标系对面积计算的影响比想象中大海岸带数据最常用的坐标系有两类一类是CGCS2000或WGS84地理坐标系以经纬度存储另一类是投影坐标系比如高斯-克吕格投影或Albers等积投影。直接在地理坐标系下计算矢量面面积单位是平方度结果会随纬度变化同一块滩涂在北方的统计值可能比南方小百分之十几这种误差在年际对比中会完全掩盖真实变化。正确做法是在读取后立刻查看坐标系统如果显示EPSG:4326或EPSG:4499这类地理坐标系后续分析必须转成等积投影。3. 用7z、Python和QGIS把RAR包转成可分析数据3.1 命令行解压先列目录、再按年定点解压在Windows下双击RAR包释放大量文件时WinRAR界面经常长时间无响应在Linux或macOS服务器上工作就更需要用命令行。7-Zip的命令行工具7z可以在不解压的情况下查看包内结构也支持按目录选择性地解压。# 列目录树确认内部结构和年度文件夹数量 7z l 18N_coastal_tidal_wetland_1989_2020.rar # 只解压2020年这一期做试读避免整包解压占用磁盘 7z e 18N_coastal_tidal_wetland_1989_2020.rar 2020/* -o./2020_test/第一条命令列出压缩包内所有文件路径第二条命令把2020年子目录下所有内容释放到单独目录里做试读。注意7z e会把所有文件平铺到一个目录如果文件名重复会发生覆盖需要保留完整路径时改用7z x# 保留目录结构完整解压适合正式使用 7z x 18N_coastal_tidal_wetland_1989_2020.rar -o./wetland_data/这里的逻辑是先看后解先把目录结构和样本文件确认一遍再全量解压减少无效解压对磁盘空间的消耗。3.2 用GeoPandas读取Shp检查几何完整性和分类字段读取矢量数据并做基本体检的Python脚本如下import geopandas as gpd fp ./wetland_data/2020/coastal_wetland_2020.shp gdf gpd.read_file(fp, encodingutf-8) print(坐标系统:, gdf.crs) print(要素数量:, len(gdf)) print(字段列表:, list(gdf.columns)) print(gdf[TYPE].value_counts()) # 检查空几何空几何会在面积计算时直接报错 invalid gdf.geometry.is_empty.sum() print(f空几何数量: {invalid})这段脚本要回答三个问题当前数据是经纬度还是投影坐标、有没有空记录、分类字段取值是否在预期范围内。gdf.crs输出为EPSG:4326说明是经纬度坐标后续面积计算必须做投影转换输出为EPSG:4499这类CGCS2000地理坐标系同样需要转换。空几何数量大于0时必须先用gdf gdf[~gdf.geometry.is_empty]过滤否则后续合并年份时会频繁报错。3.3 用Rasterio读取Tif栅格验证波段数和空间范围如果压缩包内是栅格成果用rasterio读取头信息的效率远高于直接拖进GISimport rasterio with rasterio.open(./wetland_data/1989/coastal_wetland_1989.tif) as src: print(坐标系统:, src.crs) print(空间范围:, src.bounds) print(行列数:, src.height, src.width) print(波段数:, src.count) print(像元类型:, src.dtypes)这里最需要关注的是像元类型和空间范围。分类结果通常存成uint80代表非湿地、1代表裸滩、2代表草滩、3代表红树林如果是float32则说明存的是连续型变量比如水体指数或植被指数。空间范围方面把相邻两个年份的bounds打印出来对比如果某一年和其他年份差异明显说明当年数据做了额外裁剪做逐年叠加分析时需要统一边界。3.4 QGIS里把数据落到地图上做人工复核批量处理之前在QGIS里手工加载几个年份的图层叠加天地图或Esri影像底图看细节能发现脚本难以察觉的问题。重点检查三处图层坐标系是否与底图一致、比例尺在1:50000附近时滩涂边界是否平滑、以及边界与影像上的潮沟、养殖塘轮廓是否吻合。如果加载后发现图层整体偏移几十到几百米优先怀疑坐标系基准面差异。中国海岸数据常用CGCS2000或西安80基准在线底图通常是WGS84两者在部分地区偏差几十米到百米级这不一定代表数据精度差只是坐标系转换参数没设置正确。在QGIS图层属性里把“转换设置”从默认改成精确转换多数情况下能消除这个偏移。4. 数据使用中会遇到的三个真实问题偏移、字典缺失、时序跳跃4.1 图层与底图对不上先查基准面再查配准打开数据后最常见的现象是滩涂边界与在线影像整体错开几十米方向一致、量级稳定。遇到这种情况先别怀疑数据质量大概率是基准面不同的系统偏差。查看shp元数据中记录的基准如果是GCS_Krasovsky_1940或Xian_1980就需要从该基准转换到WGS84或CGCS2000。QGIS里的处理路径是图层属性、数据源、坐标系转换设置手动指定原坐标系到目标坐标系的转换参数并勾选精确转换。Python环境直接使用pyprojfrom pyproj import Transformer # 从CGCS2000地理坐标系转到WGS84经纬度 trans Transformer.from_crs(EPSG:4499, EPSG:4326, always_xyTrue) lon, lat trans.transform(1234567.89, 4123456.78) print(lon, lat)务必保留always_xyTrue否则Transformer默认按纬度、经度顺序返回结果交换后的坐标点会落到完全错误的位置。另外要明确这类转换解决的是系统性偏移如果同一期数据在不同位置偏移方向和大小不一致有的是向岸偏、有的向海偏那是当年数据配准时引入的局部误差需要采集控制点做局部校正两种误差不能混为一谈。4.2 分类字段只有数字编号怎么恢复分类含义如果打开后看到CLASS_ID或TYPE字段里是0、1、2、4、5这样的整数而元数据文档丢失分类字典就断了。先不要猜按下面的顺序排查import geopandas as gpd gdf gpd.read_file(./wetland_data/1989/coastal_wetland_1989.shp) print(gdf[CLASS_ID].unique())输出结果里如果类别编号不连续说明分类体系里定义过但该年份没有出现的类别或者编码发生过调整。此时有两个可靠的恢复路径一是在压缩包里找附带Excel或csv格式的属性对照表二是参考国家标准《湿地分类》GB/T 24708-2009该标准把滨海湿地分成浅海水域、潮下水生层、珊瑚礁、岩石海岸、潮间沙石海滩、潮间淤泥海滩、潮间盐水沼泽、红树林等多个类型许多公开数据集直接沿用这套编号。实在无法确定时用同一年的Landsat真彩色影像叠加到图层上目视判断裸滩、植被、红树林三大类也能恢复出八成信息。4.3 某一年面积突变是真实变化还是数据质量异常时间序列上某个年份出现明显面积跳跃处理方式决定了结论是否可信。自然的围垦和淤积过程通常表现为连续几年的渐变极少发生单一年份的剧烈突变。符合这种规律的面积跳跃先怀疑影像质量问题再怀疑分类错误。import pandas as pd # area_series 是每年每个类别的面积表行为年份列为类别 area_series ( all_data.groupby([year, CLASS_ID])[area_km2].sum().unstack() ) # 只对缺失值做线性插值 area_filled area_series.interpolate(limit_directionboth)这段插值代码的适用范围很窄只能用于云量覆盖、条带缺失等随机性数据空洞不能用于插补真实发生的围填海面积减少。判断标准是看突变年份之后两年的数值是否回归到突变前水平如果回归了说明当年数据存在质量问题如果持续偏低且走势平稳那才是真实的滩涂减少。两种情况的定性要在分析报告里明确区分不能混在一起讨论。5. 做一次完整的滩涂面积变化盘点附可复用脚本5.1 批量读取全部年度矢量并统一坐标系把1989到2020年所有年度的Shapefile读进同一个GeoDataFrame是后续所有统计的前提import glob import geopandas as gpd import pandas as pd files sorted(glob.glob(./wetland_data/*/coastal_wetland_*.shp)) frames [] for fp in files: year int(fp.split(/)[-1].split(_)[2]) gdf gpd.read_file(fp) gdf[year] year frames.append(gdf) all_data gpd.GeoDataFrame(pd.concat(frames, ignore_indexTrue))这段脚本的核心是年份提取逻辑fp.split(_)[2]依赖文件名的固定格式。如果你拿到的版本命名规则不同需要调整下标或者改用正则表达式提取年份否则年份信息会整体错位。5.2 等积投影和面积统计的取舍统一转到等积投影后再计算面积。覆盖中国全境时用WGS 84/World BehrmannEPSG:6933可以避开分带问题一套参数跑完所有年份all_data all_data.to_crs(EPSG:6933) all_data[area_km2] all_data.geometry.area / 1e6 annual all_data.groupby([year, CLASS_ID])[area_km2].sum()如果项目要求高精度可以按海岸带分段改用UTM对应分带投影比如渤海用UTM 50N、东海用UTM 51N。但32年数据跨多个分带分段投影会引入分带边界上的人为断裂批量运算时不如统一使用等积投影稳定。5.3 输出年度面积表和变化率供下一步分析按区域字段对面积做聚合输出成标准CSV方便后续进入Excel或BI工具regional all_data.groupby([year, zone])[area_km2].sum() regional.unstack().to_csv(wetland_area_by_zone.csv) delta annual.unstack().diff().dropna() delta.to_csv(wetland_area_change_by_year.csv)第一张表是各区域逐年面积第二张表是相邻年份差值两张表配合使用能快速定位面积突变年份。输出的CSV里行列索引会自动带上年份和类别码直接用Excel透视表即可做进一步筛选。最后把all_data按年份切片输出当年的空间变化图斑保存成单独的GeoJSON或Shp就能与最新的围填海项目库做空间叠置比对把数据产品真正接进业务判断流程。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进