简介本资源面向从事遥感与海洋、气候变化研究的科研人员及学生聚焦Sentinel-3测高Level-2数据的读取与处理。内容围绕NetCDF格式展开涵盖利用xarray、cf-python解析海面高度、冰面高度及质量控制标志并延伸至时间序列分析、地图绘制、数据校正与结果保存等环节帮助读者掌握从原始回波产品到可用科学数据的完整流程。资源包共41个文件以py脚本、sample样例、head与master版本文件为主另含shp、shx、prj、geojson等地理空间数据及csv、description等辅助说明压缩包约2.68MB结构紧凑便于按模块查阅。目前已有1164人学习下载适合需要快速上手Sentinel-3测高数据处理、理解CF规范元数据并开展海平面与水位分析的研究者参考。1. 从一包 Sentinel-3 测高数据说起为什么 L2 才是真正能干活的层级第一次拿到 Sentinel3AltimetryL2.zip 的人多半会先愣一下里面既有 NetCDF 数据又混着DataDownload.py、code_elevation.py、ganga.shp、basin.geojson还有一个完整的.git目录。这不是一个单纯的示例数据包而是一套已经跑通过「下载—裁剪—统计—出图」链路的工程快照。Sentinel-3 是欧空局的对地观测卫星测高仪Altimeter是它的核心载荷之一能给出全球海面高度、冰盖高程、河流湖泊水位。L2 级产品是雷达回波经过初步处理后的成果海面高度、有效波高、后向散射系数、质量控制标志都在里面通常以 NetCDF 格式存储。这份资源解决的就是「怎么把 L2 从原始网格变成自己研究区里能用的水位/高程序列」——适合做海洋、水文、气候变化的人也适合想学遥感数据处理工程化的新手。2. 拆开压缩包文件清单背后的处理链路与选型逻辑2.1 每个文件在链路里的位置先把包里的东西按角色分一下不然后面写脚本会迷路。文件/目录角色说明Sentinel3AltimetryL2.zip数据主体压缩后的 L2 产品解压后是 NetCDFDataDownload.py下载脚本负责按时间/区域拉取数据code_elevation.py核心处理读取 NetCDF、提取变量、裁剪、统计FileForValues.csv中间结果从 NetCDF 抽出的表格化数值ganga.shp/.shx/.prj矢量边界研究区或流域边界basin.geojson矢量边界同上GeoJSON 格式方便 Python 读csv to shape.py格式转换CSV 转 Shapefile用于出图.git/版本痕迹说明这套脚本是迭代过的选型上作者没有用cf-python而是走xarray netCDF4路线这是常见做法xarray的接口更像 Pandas处理多维数组和坐标方便cf-python虽然能解析 CF 规范但安装重、依赖多快速迭代时反而拖节奏。矢量部分用geopandas读basin.geojson比直接读 Shapefile 少一层编码坑。2.2 环境准备把依赖一次装齐不要一个个pip install容易漏。建议先建虚拟环境再按下面这组装。xarray负责数据模型netCDF4是底层引擎geopandas处理边界cartopy出图dask在数据大时做分块。python -m venv s3env source s3env/bin/activate # Windows 用 s3env\Scripts\activate pip install xarray netCDF4 dask geopandas shapely cartopy matplotlib pandas参数说明netCDF4必须装否则xarray.open_dataset会回退到 scipy 引擎读 Sentinel-3 的组结构容易报错dask不是必须但当你要批量读几十个 L2 文件时chunks{}能让内存不爆。装完用python -c import xarray; print(xarray.__version__)验证一下版本低于 2023 的建议升级老版本对 CF 坐标的解析有玄学问题。2.3 读取 L2 并看清变量结构Sentinel-3 L2 的 NetCDF 不是平铺的变量常按组存放。先别急着取SSH先看结构。import xarray as xr # 打开单个 L2 文件group 参数按需指定 ds xr.open_dataset(Sentinel3AltimetryL2.nc, groupdata_01) # 打印所有变量和坐标确认命名 print(ds) print(ds.data_vars) print(ds.coords)逻辑说明groupdata_01是因为部分 L2 产品把测量数据放在子组里不指定会读到空壳。data_vars列出所有物理量常见的有ssh海面高度、swh有效波高、sig0后向散射、quality_flags。参数上如果文件是多个轨道段拼接的ds会带time和lat/lon坐标直接ds.ssh就能拿到 DataArray。失败时先看报错里有没有KeyError多半是组名写错用ncdump -h或xr.open_dataset(path).groups查一下真实组名。2.4 用边界裁剪从全球网格到研究区拿到全球数据后下一步是按basin.geojson裁剪。不要用ds.where硬套经纬度矩形边界不规则时会带进大量无效点。import geopandas as gpd import rioxarray # 提供 clip 能力 import xarray as xr ds xr.open_dataset(Sentinel3AltimetryL2.nc, groupdata_01) # 给 DataArray 写 CRSSentinel-3 通常是 WGS84 ds ds.rio.write_crs(EPSG:4326) basin gpd.read_file(basin.geojson) # 裁剪dropTrue 去掉边界外的坐标 clipped ds.rio.clip(basin.geometry, basin.crs, dropTrue) # 导出为 CSV对应 FileForValues.csv 的角色 df clipped[[ssh, swh, lat, lon, time]].to_dataframe().reset_index() df.to_csv(FileForValues.csv, indexFalse)逻辑说明rio.write_crs是必须的否则clip不知道数据坐标系会直接报MissingCRS。basin.crs要和数据一致不一致先用to_crs(EPSG:4326)转。to_dataframe()会把多维压成表格适合后续统计。参数上dropTrue能减少内存但如果你要保留完整网格做插值就设False。这一步的坑是time坐标可能是datetime64导出 CSV 后变成字符串读回来要pd.to_datetime转一下。3. 从 NetCDF 到 CSV 再到 Shapefile把水位序列落到地图上3.1 质量控制先过滤再统计L2 数据带quality_flags不过滤直接算均值结果会被异常点带偏。常见做法是按位掩码筛。import numpy as np # 假设 quality_flags 是整数位掩码bit 0 表示有效 valid (clipped[quality_flags] 1) 1 clean clipped.where(valid, dropTrue) # 统计有效点的平均海面高度 mean_ssh float(clean[ssh].mean().values) print(f研究区平均 SSH: {mean_ssh:.3f} m)逻辑说明 1是取最低位不同产品位定义不同先查手册确认哪一位代表「有效测量」。where(valid, dropTrue)会把无效点直接删掉而不是置 NaN后续mean不会跳过。参数上如果quality_flags是浮点先astype(int)。这一步翻车最多的是位定义搞反把无效当有效结果 SSH 出现几十米的离谱值回头查半天。3.2 CSV 转 Shapefile出图前的最后一步csv to shape.py干的就是把带经纬度的 CSV 变成点 Shapefile方便在 QGIS 或 ArcGIS 里叠加底图。import pandas as pd import geopandas as gpd from shapely.geometry import Point df pd.read_csv(FileForValues.csv) df[time] pd.to_datetime(df[time]) # 用经纬度构造点几何 geometry [Point(xy) for xy in zip(df[lon], df[lat])] gdf gpd.GeoDataFrame(df, geometrygeometry, crsEPSG:4326) # 写出 Shapefile注意字段名不要超 10 字符 gdf.to_file(ganga.shp, driverESRI Shapefile)逻辑说明Point(xy)逐个构造点数据量大时慢可以用gpd.points_from_xy(df.lon, df.lat)替代快一个量级。crsEPSG:4326必须写否则.prj文件为空别人打开会问「这数据在哪」。参数上Shapefile 字段名限制 10 个字符quality_flags这种长名会被截断导出前先rename成短名。这一步的坑是中文路径to_file遇到中文目录可能报编码错换英文路径最稳。3.3 时间序列与地图绘制有了干净的点数据就可以做时间序列和空间分布图。时间序列看某个网格点的 SSH 变化地图看整个研究区的空间格局。import matplotlib.pyplot as plt # 时间序列取最近的一个有效点 ts clean[ssh].isel(lat0, lon0).to_series() ts.plot(titleSSH Time Series) plt.savefig(ssh_timeseries.png, dpi150) # 空间分布散点图颜色映射 SSH fig, ax plt.subplots(figsize(8, 6)) sc ax.scatter(df[lon], df[lat], cdf[ssh], cmapviridis, s5) plt.colorbar(sc, labelSSH (m)) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) plt.savefig(ssh_map.png, dpi150)逻辑说明isel是按索引取适合已知网格位置如果按经纬度取用sel(lat..., lon..., methodnearest)。to_series()把 DataArray 转成带时间索引的 Seriesplot直接出图。参数上dpi150够用投稿再调 300。地图散点用s5避免点重叠成块。失败时看clean是否为空空的话说明裁剪区域和轨道没交集换时间段或换边界。4. 避坑与排查那些让我重跑一整天的细节4.1 现象open_dataset报HDF Error文件读不进去原因Sentinel-3 L2 是 NetCDF-4/HDF5 格式但下载不完整或压缩包解压时损坏文件头缺失。解决先用ncdump -h试读报错就重新下载。DataDownload.py里如果有断点续传逻辑检查.part文件是否被误当完整文件。4.2 现象裁剪后点数骤减几乎为空原因basin.geojson的坐标系不是 WGS84或者经纬度顺序写反GeoJSON 是[lon, lat]容易和[lat, lon]混。解决print(basin.crs)确认不是 4326 就to_crs再print(basin.total_bounds)看范围是否合理。4.3 现象quality_flags全为 0过滤后没数据原因部分 L2 产品的质量控制放在单独变量或属性里不在主组。解决ds.attrs和ds.data_vars都翻一遍找quality相关字段或者查产品手册确认位定义别默认 bit 0。4.4 现象CSV 转 Shapefile 后属性表乱码原因Shapefile 的 DBF 默认编码是 Latin-1中文或特殊字符会乱。解决to_file(..., encodingutf-8)或者导出前把非 ASCII 字段删掉。更稳的做法是中间用 GeoPackage最后再转 Shapefile。4.5 现象批量处理时内存爆掉原因一次open_dataset多个文件xarray默认把数据全读进内存。解决加chunks{time: 100}启用 dask 分块或者用xr.open_mfdataset时设combineby_coords让 dask 惰性加载。处理完及时ds.close()别攒着。5. 进阶把单文件脚本改成可复用的批处理与验证习惯单文件跑通只是开始真正干活时你要面对几十上百个 L2 文件。我一般会把code_elevation.py里的逻辑抽成函数用xr.open_mfdataset一次读多个再按时间分组统计。下面这个模式我用了很久稳定且好排查。import glob import xarray as xr import pandas as pd def process_batch(pattern, basin_path): files sorted(glob.glob(pattern)) # 惰性加载按时间拼接 ds xr.open_mfdataset(files, combineby_coords, chunks{time: 200}) ds ds.rio.write_crs(EPSG:4326) basin gpd.read_file(basin_path) clipped ds.rio.clip(basin.geometry, basin.crs, dropTrue) # 质量控制 valid (clipped[quality_flags] 1) 1 clean clipped.where(valid, dropTrue) # 按时间求区域均值 ts clean[ssh].mean(dim[lat, lon]).to_series() return ts ts process_batch(Sentinel3AltimetryL2_*.nc, basin.geojson) ts.to_csv(ssh_regional_timeseries.csv)逻辑说明open_mfdataset的combineby_coords会按坐标自动拼接适合时间序列文件。chunks{time: 200}控制每块大小太大内存高太小调度开销大200 是经验值。mean(dim[lat,lon])得到区域平均时间序列比逐点统计更稳。参数上如果文件时间有重叠加compatoverride避免冲突。验证方法上我习惯做三件事一是用ncdump -h对比原始文件和输出文件的变量维度确认没丢坐标二是抽一个已知浮标站的位置看 SSH 序列和浮标观测的趋势是否一致偏差在几十厘米内算正常三是把FileForValues.csv的行数和裁剪后有效点数对一下数量级不对就回头查过滤条件。这套流程走下来基本能保证结果可复现。从那以后我每次处理新的 L2 数据都强制先跑一遍ncdump -h和边界范围打印再动裁剪和统计。这个习惯帮我省了至少三次重跑一整天的时间。希望帮到你。本文还有配套的精品资源点击获取 Sentinel-3测高L2数据处理：从NetCDF到Shapefile的工程化实战
简介本资源面向从事遥感与海洋、气候变化研究的科研人员及学生聚焦Sentinel-3测高Level-2数据的读取与处理。内容围绕NetCDF格式展开涵盖利用xarray、cf-python解析海面高度、冰面高度及质量控制标志并延伸至时间序列分析、地图绘制、数据校正与结果保存等环节帮助读者掌握从原始回波产品到可用科学数据的完整流程。资源包共41个文件以py脚本、sample样例、head与master版本文件为主另含shp、shx、prj、geojson等地理空间数据及csv、description等辅助说明压缩包约2.68MB结构紧凑便于按模块查阅。目前已有1164人学习下载适合需要快速上手Sentinel-3测高数据处理、理解CF规范元数据并开展海平面与水位分析的研究者参考。1. 从一包 Sentinel-3 测高数据说起为什么 L2 才是真正能干活的层级第一次拿到 Sentinel3AltimetryL2.zip 的人多半会先愣一下里面既有 NetCDF 数据又混着DataDownload.py、code_elevation.py、ganga.shp、basin.geojson还有一个完整的.git目录。这不是一个单纯的示例数据包而是一套已经跑通过「下载—裁剪—统计—出图」链路的工程快照。Sentinel-3 是欧空局的对地观测卫星测高仪Altimeter是它的核心载荷之一能给出全球海面高度、冰盖高程、河流湖泊水位。L2 级产品是雷达回波经过初步处理后的成果海面高度、有效波高、后向散射系数、质量控制标志都在里面通常以 NetCDF 格式存储。这份资源解决的就是「怎么把 L2 从原始网格变成自己研究区里能用的水位/高程序列」——适合做海洋、水文、气候变化的人也适合想学遥感数据处理工程化的新手。2. 拆开压缩包文件清单背后的处理链路与选型逻辑2.1 每个文件在链路里的位置先把包里的东西按角色分一下不然后面写脚本会迷路。文件/目录角色说明Sentinel3AltimetryL2.zip数据主体压缩后的 L2 产品解压后是 NetCDFDataDownload.py下载脚本负责按时间/区域拉取数据code_elevation.py核心处理读取 NetCDF、提取变量、裁剪、统计FileForValues.csv中间结果从 NetCDF 抽出的表格化数值ganga.shp/.shx/.prj矢量边界研究区或流域边界basin.geojson矢量边界同上GeoJSON 格式方便 Python 读csv to shape.py格式转换CSV 转 Shapefile用于出图.git/版本痕迹说明这套脚本是迭代过的选型上作者没有用cf-python而是走xarray netCDF4路线这是常见做法xarray的接口更像 Pandas处理多维数组和坐标方便cf-python虽然能解析 CF 规范但安装重、依赖多快速迭代时反而拖节奏。矢量部分用geopandas读basin.geojson比直接读 Shapefile 少一层编码坑。2.2 环境准备把依赖一次装齐不要一个个pip install容易漏。建议先建虚拟环境再按下面这组装。xarray负责数据模型netCDF4是底层引擎geopandas处理边界cartopy出图dask在数据大时做分块。python -m venv s3env source s3env/bin/activate # Windows 用 s3env\Scripts\activate pip install xarray netCDF4 dask geopandas shapely cartopy matplotlib pandas参数说明netCDF4必须装否则xarray.open_dataset会回退到 scipy 引擎读 Sentinel-3 的组结构容易报错dask不是必须但当你要批量读几十个 L2 文件时chunks{}能让内存不爆。装完用python -c import xarray; print(xarray.__version__)验证一下版本低于 2023 的建议升级老版本对 CF 坐标的解析有玄学问题。2.3 读取 L2 并看清变量结构Sentinel-3 L2 的 NetCDF 不是平铺的变量常按组存放。先别急着取SSH先看结构。import xarray as xr # 打开单个 L2 文件group 参数按需指定 ds xr.open_dataset(Sentinel3AltimetryL2.nc, groupdata_01) # 打印所有变量和坐标确认命名 print(ds) print(ds.data_vars) print(ds.coords)逻辑说明groupdata_01是因为部分 L2 产品把测量数据放在子组里不指定会读到空壳。data_vars列出所有物理量常见的有ssh海面高度、swh有效波高、sig0后向散射、quality_flags。参数上如果文件是多个轨道段拼接的ds会带time和lat/lon坐标直接ds.ssh就能拿到 DataArray。失败时先看报错里有没有KeyError多半是组名写错用ncdump -h或xr.open_dataset(path).groups查一下真实组名。2.4 用边界裁剪从全球网格到研究区拿到全球数据后下一步是按basin.geojson裁剪。不要用ds.where硬套经纬度矩形边界不规则时会带进大量无效点。import geopandas as gpd import rioxarray # 提供 clip 能力 import xarray as xr ds xr.open_dataset(Sentinel3AltimetryL2.nc, groupdata_01) # 给 DataArray 写 CRSSentinel-3 通常是 WGS84 ds ds.rio.write_crs(EPSG:4326) basin gpd.read_file(basin.geojson) # 裁剪dropTrue 去掉边界外的坐标 clipped ds.rio.clip(basin.geometry, basin.crs, dropTrue) # 导出为 CSV对应 FileForValues.csv 的角色 df clipped[[ssh, swh, lat, lon, time]].to_dataframe().reset_index() df.to_csv(FileForValues.csv, indexFalse)逻辑说明rio.write_crs是必须的否则clip不知道数据坐标系会直接报MissingCRS。basin.crs要和数据一致不一致先用to_crs(EPSG:4326)转。to_dataframe()会把多维压成表格适合后续统计。参数上dropTrue能减少内存但如果你要保留完整网格做插值就设False。这一步的坑是time坐标可能是datetime64导出 CSV 后变成字符串读回来要pd.to_datetime转一下。3. 从 NetCDF 到 CSV 再到 Shapefile把水位序列落到地图上3.1 质量控制先过滤再统计L2 数据带quality_flags不过滤直接算均值结果会被异常点带偏。常见做法是按位掩码筛。import numpy as np # 假设 quality_flags 是整数位掩码bit 0 表示有效 valid (clipped[quality_flags] 1) 1 clean clipped.where(valid, dropTrue) # 统计有效点的平均海面高度 mean_ssh float(clean[ssh].mean().values) print(f研究区平均 SSH: {mean_ssh:.3f} m)逻辑说明 1是取最低位不同产品位定义不同先查手册确认哪一位代表「有效测量」。where(valid, dropTrue)会把无效点直接删掉而不是置 NaN后续mean不会跳过。参数上如果quality_flags是浮点先astype(int)。这一步翻车最多的是位定义搞反把无效当有效结果 SSH 出现几十米的离谱值回头查半天。3.2 CSV 转 Shapefile出图前的最后一步csv to shape.py干的就是把带经纬度的 CSV 变成点 Shapefile方便在 QGIS 或 ArcGIS 里叠加底图。import pandas as pd import geopandas as gpd from shapely.geometry import Point df pd.read_csv(FileForValues.csv) df[time] pd.to_datetime(df[time]) # 用经纬度构造点几何 geometry [Point(xy) for xy in zip(df[lon], df[lat])] gdf gpd.GeoDataFrame(df, geometrygeometry, crsEPSG:4326) # 写出 Shapefile注意字段名不要超 10 字符 gdf.to_file(ganga.shp, driverESRI Shapefile)逻辑说明Point(xy)逐个构造点数据量大时慢可以用gpd.points_from_xy(df.lon, df.lat)替代快一个量级。crsEPSG:4326必须写否则.prj文件为空别人打开会问「这数据在哪」。参数上Shapefile 字段名限制 10 个字符quality_flags这种长名会被截断导出前先rename成短名。这一步的坑是中文路径to_file遇到中文目录可能报编码错换英文路径最稳。3.3 时间序列与地图绘制有了干净的点数据就可以做时间序列和空间分布图。时间序列看某个网格点的 SSH 变化地图看整个研究区的空间格局。import matplotlib.pyplot as plt # 时间序列取最近的一个有效点 ts clean[ssh].isel(lat0, lon0).to_series() ts.plot(titleSSH Time Series) plt.savefig(ssh_timeseries.png, dpi150) # 空间分布散点图颜色映射 SSH fig, ax plt.subplots(figsize(8, 6)) sc ax.scatter(df[lon], df[lat], cdf[ssh], cmapviridis, s5) plt.colorbar(sc, labelSSH (m)) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) plt.savefig(ssh_map.png, dpi150)逻辑说明isel是按索引取适合已知网格位置如果按经纬度取用sel(lat..., lon..., methodnearest)。to_series()把 DataArray 转成带时间索引的 Seriesplot直接出图。参数上dpi150够用投稿再调 300。地图散点用s5避免点重叠成块。失败时看clean是否为空空的话说明裁剪区域和轨道没交集换时间段或换边界。4. 避坑与排查那些让我重跑一整天的细节4.1 现象open_dataset报HDF Error文件读不进去原因Sentinel-3 L2 是 NetCDF-4/HDF5 格式但下载不完整或压缩包解压时损坏文件头缺失。解决先用ncdump -h试读报错就重新下载。DataDownload.py里如果有断点续传逻辑检查.part文件是否被误当完整文件。4.2 现象裁剪后点数骤减几乎为空原因basin.geojson的坐标系不是 WGS84或者经纬度顺序写反GeoJSON 是[lon, lat]容易和[lat, lon]混。解决print(basin.crs)确认不是 4326 就to_crs再print(basin.total_bounds)看范围是否合理。4.3 现象quality_flags全为 0过滤后没数据原因部分 L2 产品的质量控制放在单独变量或属性里不在主组。解决ds.attrs和ds.data_vars都翻一遍找quality相关字段或者查产品手册确认位定义别默认 bit 0。4.4 现象CSV 转 Shapefile 后属性表乱码原因Shapefile 的 DBF 默认编码是 Latin-1中文或特殊字符会乱。解决to_file(..., encodingutf-8)或者导出前把非 ASCII 字段删掉。更稳的做法是中间用 GeoPackage最后再转 Shapefile。4.5 现象批量处理时内存爆掉原因一次open_dataset多个文件xarray默认把数据全读进内存。解决加chunks{time: 100}启用 dask 分块或者用xr.open_mfdataset时设combineby_coords让 dask 惰性加载。处理完及时ds.close()别攒着。5. 进阶把单文件脚本改成可复用的批处理与验证习惯单文件跑通只是开始真正干活时你要面对几十上百个 L2 文件。我一般会把code_elevation.py里的逻辑抽成函数用xr.open_mfdataset一次读多个再按时间分组统计。下面这个模式我用了很久稳定且好排查。import glob import xarray as xr import pandas as pd def process_batch(pattern, basin_path): files sorted(glob.glob(pattern)) # 惰性加载按时间拼接 ds xr.open_mfdataset(files, combineby_coords, chunks{time: 200}) ds ds.rio.write_crs(EPSG:4326) basin gpd.read_file(basin_path) clipped ds.rio.clip(basin.geometry, basin.crs, dropTrue) # 质量控制 valid (clipped[quality_flags] 1) 1 clean clipped.where(valid, dropTrue) # 按时间求区域均值 ts clean[ssh].mean(dim[lat, lon]).to_series() return ts ts process_batch(Sentinel3AltimetryL2_*.nc, basin.geojson) ts.to_csv(ssh_regional_timeseries.csv)逻辑说明open_mfdataset的combineby_coords会按坐标自动拼接适合时间序列文件。chunks{time: 200}控制每块大小太大内存高太小调度开销大200 是经验值。mean(dim[lat,lon])得到区域平均时间序列比逐点统计更稳。参数上如果文件时间有重叠加compatoverride避免冲突。验证方法上我习惯做三件事一是用ncdump -h对比原始文件和输出文件的变量维度确认没丢坐标二是抽一个已知浮标站的位置看 SSH 序列和浮标观测的趋势是否一致偏差在几十厘米内算正常三是把FileForValues.csv的行数和裁剪后有效点数对一下数量级不对就回头查过滤条件。这套流程走下来基本能保证结果可复现。从那以后我每次处理新的 L2 数据都强制先跑一遍ncdump -h和边界范围打印再动裁剪和统计。这个习惯帮我省了至少三次重跑一整天的时间。希望帮到你。本文还有配套的精品资源点击获取