ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

NSCAT L3海面风场数据处理全指南:从NetCDF读取到绘图

NSCAT L3海面风场数据处理全指南:从NetCDF读取到绘图 手上刚好整理过一批 NSCAT 3 级风场数据翻了翻当时的处理脚本和踩坑记录干脆把这些东西梳理成一篇完整的说明。无论你是刚接触散射计数据、准备做海面风场分析还是想在论文里引用这套历史数据这篇内容应该能帮你少走不少弯路。1. 这套数据是什么来头NSCAT 任务背景与数据定位NSCATNASA Scatterometer美国宇航局散射计搭载在1996年8月发射的日本先进地球观测卫星ADEOS-I上是当时全球第一颗业务化运行的Ku波段散射计。它通过测量海面粗糙度引起的雷达后向散射系数变化反演海面10米高度处的风矢量覆盖全球无冰海洋面积约78%。1997年6月ADEOS-I卫星因太阳能帆板故障导致整星失效NSCAT实际有效观测时间只有大约9个多月。虽然寿命短但它为后续QuikSCAT、ASCAT等散射计提供了关键的技术验证也是1990年代中期唯一能够提供高分辨率全球海洋风场的星载传感器。JPL喷气推进实验室基于NSCAT L2B轨道风矢量产品制作了这套3级每日网格化浏览图像。它本质上是一种快速可视化产品把不规则分布的轨道风场重采样到等经纬度网格生成每日全球风场分布图输出为图像和NetCDF数据两种形式。因为NSCAT的观测窗口比较特殊而且L2B数据使用门槛较高很多利用NSCAT研究热带气旋、海洋环流、海气相互作用的论文都会直接采用这套L3产品。你需要明确一个关键点NSCAT L3浏览图像的定位是快速查看和统计分析基础。它的空间分辨率是1°×1°时间分辨率是每日一双升轨和降轨分别合成适合分析天气尺度以上的风场特征不适合做中尺度锋面或台风内部精细结构的个例研究。后面我会详细解释为什么网格化之后分辨率会变粗以及哪些应用场景下你需要回到L2B。数据覆盖时间从1996年9月15日到1997年6月29日每日一个文件全球范围覆盖。下载途径一般是NASA物理海洋学分布式数据档案中心PO.DAAC现在通过Earthdata Search界面检索NS CAT Level 3即可找到。老用户在1998-2005年间常用的FTP路径已经发生了变化HTTPS和OPeNDAP接口现在是主要的数据访问方式。2. 文件里到底有什么NetCDF结构、变量含义与元数据解读这套3级产品最常用的格式是NetCDF文件命名大致遵循s-nscat-l3-daily-YYYYMMDD-...nc的形式。打开文件后你首先会看到维度定义和几个核心变量。对于不熟悉NetCDF的读者可以把它理解成一种自描述的数组容器变量名、单位、缺失值都写在文件内部比纯二进制和ASCII要友好得多。核心变量大致如下变量名维度单位含义wind_speed(lat, lon)m/s海面10米高度风速wind_direction(lat, lon)degree风的来向气象学角度u_wind/v_wind(lat, lon)m/s风矢量的纬向和经向分量n_measurements(lat, lon)count参与平均的观测次数time1days since 1996-01-01参考时间latitude/longitude1Ddegree网格中心坐标quality_flag(lat, lon)-质量标记rain_flag(lat, lon)-降雨影响标记其中u_wind和v_wind是最便于直接使用的变量因为风场分析经常涉及散度、涡度、通量计算直接用分量比风速风向转换方便得多。wind_direction是气象学中的来向定义0度表示北风风从北往南吹90度表示东风。做风矢量图的时候要注意在数学坐标系中绘制时需要做u-wind_speed*sin(dir*pi/180)或者直接用已有的u_wind、v_wind不要自己拿速度和方向去转否则方向和速度对应不上很容易出错。n_measurements这个变量很多人会忽略但它非常有用。它表示某网格内有多少个有效观测值参与了平均数值越高该格点的风场可信度越高。如果某个网格的n_measurements等于0说明没有升轨或降轨数据覆盖。需要注意的是NSCAT的一个轨道覆盖带宽度约600公里相邻轨道之间在中低纬度存在缝隙所以不是每天每个格点都有数据。在做逐日图时这些n_measurements0的区域要设置成NaN不能当作缺测风场处理更不能把0当成风速值来画图。元数据部分记录了数据生成的版本号、输入L2B产品版本、网格化算法简介、以及JPL数据处理团队的联系方式。早期版本如V2、V3与后续重处理版本如V4之间风速和风向存在细微的系统性差异。建议你在论文方法部分明确写出你使用的是哪个版本因为不同版本的风速均值差异最多可达0.3-0.5 m/s虽然看起来小但在气候统计分析里足以改变趋势和显著性结论。文件内部还包含了一些辅助质量信息包括每个网格使用的观测时间窗、平均后向散射系数等但这些字段在不同发布版本中差异较大需要查阅当时的数据说明文档README来确定具体含义。早期JPL发布的数据说明文档往往以PDF形式附带在数据目录中比NetCDF全局属性更详细。3. 读取与预处理从NetCDF到可分析数据的关键步骤拿到NetCDF文件后最直接的处理方式是使用Python的xarray库。不仅因为xarray对NetCDF的维度、变量、坐标处理非常优雅还因为它内置了resample和groupby方法便于做多日合成和气候态计算。下面是一个典型的读取和预处理流程也是我实际处理这套数据时的标准操作。3.1 读取并生成掩码import xarray as xr import numpy as np ds xr.open_dataset(s-nscat-l3-daily-19961001-...nc) lat ds[latitude].values lon ds[longitude].values # 读取U/V分量 u ds[u_wind].squeeze() v ds[v_wind].squeeze() speed np.sqrt(u**2 v**2) # 质量与降雨标记 qual ds[quality_flag].squeeze() rain ds[rain_flag].squeeze() nmeas ds[n_measurements].squeeze() # 有效数据掩码质量标记合格、无降雨、观测数大于0 mask (qual 0) (rain 0) (nmeas 0) u u.where(mask) v v.where(mask) speed speed.where(mask)这里有几个重要的细节需要说明。第一quality_flag0通常表示质量合格非0值表示数据受各种因素污染包括陆地污染、海冰误判、极端风速不可信等。如果你忽略标记直接使用所有数据沿海岸线和高纬度地区会出现很多高风速的奇异点。第二rain_flag对于Ku波段散射计非常重要。Ku波段14 GHz信号对雨滴敏感降雨会显著改变海面粗糙度导致反演风速严重偏高尤其是大雨条件下风速可能被高估5-10 m/s。rain_flag标记了受降雨影响的观测在气候态统计分析中应该剔除。第三nmeas0这个条件看上去简单却是最容易翻车的地方。有些处理工具在网格化填充时会把空值填为0如果你不设这个条件直接把所有风速为0的区域当成静风在全球图上就会看到轨道缝隙处出现大片0风速带这完全不符合实际海况。3.2 建立陆地/海冰掩码虽然NSCAT L3产品已经标注了大部分陆地像元但近岸网格受陆地回波污染的情况仍会出现。一种常见做法是使用NSCAT L3自带的land_mask或sea_ice_flag变量如果有另一种是使用外部海岸线数据生成掩码。# 如果文件内没有掩码变量可以用外部海岸线生成 # 这里使用Natural Earth的1:110m海岸线示意 import geopandas as gpd from shapely.geometry import Point world gpd.read_file(ne_110m_land.shp) lon_2d, lat_2d np.meshgrid(lon, lat) mask_land np.zeros(lon_2d.shape, dtypebool) for i in range(lat_2d.shape[0]): for j in range(lat_2d.shape[1]): pt Point(lon_2d[i,j], lat_2d[i,j]) mask_land[i,j] world.contains(pt)不过这种双重循环在1°网格上跑全球180×360倒还好但如果以后处理0.25°数据就会慢得让人崩溃。更好的做法是用regionmask库或者用cartopy的feature接口先生成岸线栅格。例如import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.io import shapereader # 用cartopy在1°网格上生成陆地掩码 import xarray as xr land_110 cfeature.NaturalEarthFeature( physical, land, 110m, edgecolorface, facecolornone) lon_b, lat_b np.meshgrid(lon, lat) land_mask np.zeros(lon_b.shape, dtypebool) # 这里示意用经纬度范围初筛实际可用shapely的contains # regionmask库一行即可完成推荐 import regionmask mask_region regionmask.defined_regions.natural_earth_v5_0_0.land_110() land_mask mask_region.mask(lon_2d, lat_2d).values land_mask ~np.isnan(land_mask) # True表示陆地我的经验是对于NSCAT这种1°浏览产品用文件自带的land_flag如果有加上一个简单的近岸缓冲比如距离海岸线0.5°以内的格点全部剔除就够用了不需要精确到像元的陆地掩码。原因是L3本身的网格化分辨率有限海岸附近的风场反演误差远大于分辨率损失你在海岸线附近看到的强风很可能是陆地回波泄漏导致的假信号。3.3 升轨与降轨数据的区分使用NSCAT的每日L3文件通常包含升轨ascending和降轨descending的分别合成变量或者通过orbit_flag标记不同轨道方向。为什么需要区分因为升轨和降轨的观测地方时不同。对NSCAT来说ADEOS卫星的轨道设计使升轨大约在地方时10:30左右经过赤道降轨大约在22:30左右。海面风场具有明显的日变化尤其是海陆风影响区域和热带对流活动区如果把两个时刻的数据混在一起做日均会引入至少1 m/s量级的内噪声。在1996-1997年的NSCAT L3版本中某些时段文件里升轨和降轨变量是分开命名的比如u_wind_asc、u_wind_desc。在做日平均风场时我通常先分别检查两个方向风速的差异如果差异超过2 m/s的区域比例大于10%说明当天天气系统活跃日均值需要谨慎解读。实际做气候态时建议升轨和降轨分别计算平均值再对两个平均场做算术平均这样等效于对日变化做了最简单的订正。4. 画风场图的实际操作一张可用的风矢量图是怎样生成的对于海洋风场来说一张好的浏览图像应当同时表达风速大小和方向。JPL原始浏览图像一般以颜色填充表示风速叠加箭头表示风矢量方向。下面以Python的matplotlib和cartopy为例给出一个可以复制使用的绘图流程并说明几个容易忽略的细节。4.1 风场填色图的绘制逻辑import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 以某一天数据为例示意假设已经从xarray中取出u、v、speed fig plt.figure(figsize(14, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.set_global() # 绘制风速填色 cf ax.pcolormesh(lon, lat, speed, cmapSpectral_r, vmin0, vmax20, shadingauto) # 叠加海岸线 ax.add_feature(cfeature.LAND, facecolorlightgray) ax.add_feature(cfeature.COASTLINE, linewidth0.5) # 叠加风矢量箭头直接绘制U/V分量 # 为了清晰每3个格点画一个箭头 q ax.quiver(lon[::3], lat[::3], u[::3, ::3], v[::3, ::3], scale400, width0.002, alpha0.7, transformccrs.PlateCarree()) cbar fig.colorbar(cf, axax, orientationhorizontal, pad0.05, shrink0.8) cbar.set_label(Wind speed (m/s)) ax.set_title(NSCAT L3 Daily Wind Vector 1996-10-01) plt.show()这段代码基本复现了JPL浏览图像的核心内容。使用Spectral_r的原因是这个colormap在0-20 m/s范围内的视觉层次分明低风速偏蓝、高风速偏红符合海洋气象学社区的习惯。如果你需要更接近JPL原始图像的配色可以用ncl_default或viridis以外的一些气象专用colormap比如cmcolor或cmocean的balance。4.2 出图时容易被忽略的坐标陷阱在绘制这套L3数据时有一个非常常见的坑文件里的longitude坐标可能是0到360也可能是-180到180不同版本不统一。如果你不做诊断直接画从0到360的经度范围在PlateCarree投影上会导致地图中心位置偏移图形看起来像被拉断了。建议读过文件后先检查lon.min()和lon.max()如果大于180利用xarray的assign_coords方法转换为-180到180或者直接对lon进行((lon 180) % 360) - 180变换。另一个坑是关于quiver的投影参数。在cartopy中当你用ax.quiver绘制U/V分量时变换需要设置transformccrs.PlateCarree()否则风矢量在非等经纬度投影下会指向错误的方位。NSCAT本身是全球覆盖的我最常用的投影是PlateCarree适合快速查看全球图和Orthographic适合极地视角。如果做区域图比如热带太平洋或北大西洋建议用Mercator或LambertConformal这时风矢量方向会随纬度变化必须依赖cartopy正确处理矢量旋转不要让箭头飞出图框。4.3 多久采样一个箭头合适很多初学者一上来就把每个格点的风矢量都画上去结果整个图形变成一片密密麻麻的黑色箭头什么都看不清。1°网格的全球图有180×360个格点每个格点上都画一个箭头显然不现实。我的经验是全球图每4-5个格点取一个箭头区域图每2个格点取一个箭头。取样的方式不要用np.arange简单隔点选取而要确认取样后仍然覆盖了关键系统比如气旋中心、急流带必要时可以先用scipy.ndimage.uniform_filter对u/v做轻微平滑再降采样。这样能保证箭头方向不是噪声。对NSCAT L3来说原始数据本身已经是日均产品网格内的风场已经被平滑过再做平滑对趋势分析影响不大但要注意不要对wind_speed填色图做过多的平滑因为浏览图像的一个核心功能是呈现空间变率如果你为了美观把风速场平滑到看不出锋面结构那就失去意义了。5. 版本坑与质量标记为什么同一区域数值和别人对不上使用NSCAT L3数据最让人头疼的问题是你得到的数值可能和别人论文里的数值存在系统偏差。这不一定是操作错误更可能是使用了不同的数据版本和不同的质量筛选策略。整理一下常见的原因方便你排查。5.1 网格化算法差异NSCAT L3历史上经历过多次重处理。最早的V1版本把轨道数据简单平均到网格V2引入了距离加权插值V3加入了海冰和降雨标记V4修正了风和海况状态的边界处理。不同版本的同一日平均风速在赤道辐合带和高纬度西风带可以相差0.5-1 m/s。如果论文里没有写明确版本号很难直接对比。JPL在1998年发布的第一版官方文档中把L3产品描述为为用户提供快速浏览能力这让很多人误以为它是研究级产品。但从后续评估来看L3在热带地区的风速偏差和均方根误差确实比同时段的L2B轨道产品要大。原因是网格化过程中的平滑效应会损失部分风速极值和梯度信息而风场反演里最需要关心的恰恰是热带气旋、冷空气爆发这类极端事件。对于气候平均意义上的风速研究L3完全够用但对于台风个例和强对流系统的风场结构分析还是应该回到L2B。5.2 降雨标记是推荐使用还是必须使用我在前文已经提到rain_flag的重要性这里再展开说说。Ku波段散射计发展早期雨的影响曾经是一个令人头疼的问题。雨滴对雷达波的散射和衰减同时起作用导致σ0显著升高。NSCAT的降雨标记是根据辐射计和散射计联合估计得到的它不仅标记了大雨也标记了中等降雨影响。数据文档中建议所有应用都使用降雨标记来排除数据但实际操作中很多人为了保留更多样本会直接忽略这个标记。我的建议是做气候态统计时必须剔除rain_flagTrue的像元在个例分析中如果要保留至少要在图上标出降雨影响区域并合理评估这部分数据的可靠性。快速估算显示热带地区的降雨影响像元比例可达10-20%如果你不剔除全球平均风速会被抬高约0.3-0.4 m/s。这个偏量看似小但在研究年际变化和趋势时足以导致结论反转。5.3 质量标记的不同位域含义NSCAT L3的quality_flag是一个位掩码字段不同数值组合对应不同的质量问题。有些版本用0表示好、1表示坏有些版本则定义了一系列位。如果你只判断数值是否等于0有可能会把可用但质量较差的像元全部剔除也可能保留了看起来很好但实际是陆地污染的像元。最稳妥的方法是在数据集中找到质量标记的定义说明通常在全局属性flag_meanings中用位运算检查每一个位# 假设某版本flag第0位表示陆地污染第1位表示海冰第2位表示降雨 bad_land (qual 1) ! 0 bad_ice (qual 2) ! 0 bad_rain (qual 4) ! 0对于NSCAT L3的多数版本直接判断qual 0是可行的简化操作但在正式分析中建议解析每一位的质量信息至少把陆地污染和海冰误判单独提取出来看看它们在空间上是否集中在某些区域。这样可以避免把系统误差当成信号。6. 从浏览图像到实际研究NSCAT L3的典型应用与统计预处理不管你是做海洋环流、海气通量、还是极地科学拿到NSCAT L3以后最常做的操作首先是验证数据合理性然后按月或按季合成气候态最后再做具体分析。这里分享两个层面的经验。6.1 数据合理性验证在信任任何数据集之前先做基本的合理性检查。对L3风场我会依次检查这几项全球风速平均值应大致在6-8 m/s之间取决于季节和是否剔除降雨区域如果全球日均风速低于5 m/s或者高于10 m/s说明可能没有正确过滤陆地/冰面像元或者大雨污染极其严重。画一张任意一天的全球风速图确认低纬度存在明显的信风带结构中纬度西风带风速较强极地东风带弱且破碎。检查北极和南极海冰边缘区的风速那里经常出现虚假的高风速海冰表面回波和海面不同散射计反演算法在海冰区域失效。用MODIS或NSIDC海冰密集度数据做掩膜能有效剥离这部分虚假信号。和浮标数据做点对点比较时要注意时空匹配窗口。NSCAT的观测瞬时值和浮标小时平均在风速上存在约1 m/s的差异风向差异约20度尤其是在低风速条件下。散射计反演在4-6 m/s以下的低风速段风向误差很大因为此时海面粗糙度主要由风浪而不是涌浪决定方向信号弱。6.2 月平均合成与缺测补插NSCAT数据只有9个多月跨了1996年9月到1997年6月四季还算齐全但没有覆盖完整年度周期。要做准气候态可以把同一月份的若干天数据做平均比如1996年10月1-31日和1996年10月的日平均再平均。但因为轨道空隙的存在某些格点一个月内可能只有一半的天数有有效观测直接平均会造成偏差。建议在做逐月平均前先统计每个格点的有效天数有效天数少于10天的格点直接标记为缺测。缺测区域的补插是另一个大坑。全球L3网格上的空值带轨道缝隙是有规律分布的中低纬度地区每天都存在。如果你对全球风场做EOF分析或计算散度这些空值会导致结果出现条带状虚假信号。常见做法是使用最优插值OI或克里金插值但我不建议在逐日尺度上做复杂的空间补插。原因很简单NSCAT数据的轨道缝隙内可能确实存在中小尺度天气系统而你用周边插值补出来的风场本质上是一种平滑猜测在后续物理量计算中可能放大噪声。如果你的目标是一个月的平均风场更好的做法是先把每天的U/V网格数据放到一个三维数组中对每个格点的时间序列求平均同时保留有效天数。不同区域的有效天数差异很大极地轨道重叠多有效天数多赤道附近轨道相邻间距大有效天数少直接用平均场绘图的可靠性在不同纬度带差别很大。出图时把有效天数也画出来能让读者对可信度有直观感受。6.3 NSCAT与后续散射计数据的时间衔接NSCAT失效后直到1999年QuikSCAT发射中间约两年的全球散射计风场存在空缺。如果你想做更长时序的风场变化分析必然要把NSCAT和QuikSCAT衔接起来。这里有个明显的系统偏差问题NSCAT和QuikSCAT虽然都是Ku波段散射计但天线配置、入射角范围、反演算法版本不同两者在相同海域的风速存在约0.3-0.5 m/s的偏差风向偏差小于5度。直接拼接使用会导致1996-1997年和1999-2009年的风速序列出现跳变。解决方法是使用重叠期数据做统计校准。NSCAT和QuikSCAT实际没有重叠观测时段因此无法直接做同平台交叉定标。常用的校准桥梁是NDBC浮标和TAO/TRITON浮标阵列。你可以在两个时间段内分别把卫星数据和浮标数据做回归然后通过浮标这个公共参考把NSCAT调整到QuikSCAT的基准。这个两步校正的过程会累积一定误差但比直接拼接要可信得多。此外NSCAT的一个独特价值在于它能同时提供VV和HH两种极化的背向散射系数。这在后续的QuikSCAT只有VV极化和ASCATVV极化为主上是没有的。不过L3浏览产品并不直接包含极化信息如果你需要研究极化比模型应该直接使用NSCAT L1.7或L2A级数据。在L3尺度上NSCAT更多被用来提供海面风场的大尺度分布特征。7. 从产品文档到实测效果的几个结论说回到JPL原始浏览图像本身。JPL在发布L3浏览图像时以彩色图形式直观展示每日全球风场。这种图像适合新闻稿、科普和快速检查但正式研究不能把JPG图片作为数据源。原因是JPEG压缩会损失色标精度每一个色标对应的风速区间在0.5-1 m/s量级同时图像本身不携带精确的空间坐标和元数据。正确做法是直接下载NetCDF数据自己在本地绘制。另一个容易被忽视的问题是PODAAC网站上的NSCAT L3产品在经过系统迁移后部分旧版文件的命名和变量结构发生了细微变化。建议下载文件后立刻用ncdump -h查看变量列表确认你要用的变量确实存在。例如早期版本的wind_speed在文件里可能叫ws后来统一改成wind_speed。如果直接套用旧脚本很容易报KeyError。最后想提醒一点NSCAT这套浏览图像虽然只有9个月数据但它捕捉到了1997-1998年厄尔尼诺事件发展初期热带太平洋风场异常的演变过程。结合1997年后期的浮标和再分析资料NSCAT风场数据仍然是研究这次强厄尔尼诺事件中大气海洋响应不可多得的观测资料。如果你正在研究那个时段的海气过程花点时间把这套L3数据处理干净往往能得到一些再分析资料里被平滑掉的真实信号。
返回列表