satellite python matplotlib

NCEP FNL数据读取分析讲解

Posted by Atmospheric Chemistry & Satellite Remote Sensing on May 23, 2023

前言

FNL 文件名的规律很简单：如：fnl_20210720_00_00.grib2，就是 2021 年 7 月 20 日 0 时的 FNL 再分析数据。

变量与垂直维度介绍

从下载数据的网站还可以了解到，FNL 再分析数据空间分辨率为 1°X1°，时间分辨率为 6 小时，其中包含温度、湿度、风速、降水等多个变量。变量说明表地址： https://rda.ucar.edu/datasets/ds083.2/#metadata/grib2.html?_do=y

高度/垂直维度的介绍，则在 https://rda.ucar.edu/datasets/ds083.2/#metadata/grib2_levels.html?_do=y

记不住也没关系，能记住或者查询到你常用的就好。待会我们跟随两个案例实际感受下怎么用。

二、读取数据的基本流程

通常来讲，读取某个数据，大致分为三步：

  1. 环境准备阶段：确定操作环境和操作语言，安装必要的库。
  2. 下载数据文件：了解文件的信息。如文件名、文件属性、数据维度、数据变量名等。
  3. 确定读取文件的参数：如数据的存放路径、数据文件的名称、需要读取变量的名称及维度等。📌 尤其需要注意大小写、中英文字符、特殊符号是否需要转义字符等。

在气象/气候领域，使用较多的数据格式有 nc、grib、hdf5、csv 等等，常用的数据处理工具有 python、NCL、R、Fortran 等等，常用的操作系统主要有 Windows 和 Linux，在需要处理大规模数据或有大规模运算需求的，常在 Linux 操作系统。

## 打印镜像里的包版本
import xarray as xr
import Nio

print("xarry version: " + xr.__version__)
print("Nio version: " + Nio.__version__)

三、两个案例，掌握读取单个 grib 格式 FNL 数据

我们以世界时 2021 年 7 月 20 日 00 时的 FNL 再分析数据为例，演示如何下载与读取，并查阅它的详细信息，最后挑选若干感兴趣的变量。

1. 新建 notebook，选择好环境

在本例中，使用 python 的 xarray 及 pyNIO 库2。

# 导入必要的库，此处仅需要xarray和pynio
import xarray as xr
import Nio

2. 下载指定的 FNL 文件，挂载进项目

本部分略过。

3. 读取文件，查看详细信息

在导入必要的库后，加入以下语句，将数据文件读取到 ds 中：

In [3]:

ds = xr.open_dataset("/home/mw/input/fnl_test2986/fnl_20210720_00_00.grib2", engine="pynio") # 读取2021年7月20日0时FNL再分析数据
ds # 查看数据

Out[3]:

👆 代码输出中，你可以看到维度信息（Dimensions）、坐标信息（Coordinates）以及变量信息（Data variables），它们可以帮助你初步了解该文件内数据的情况。

如：

  • lat_0 这个维度是从 90.0° 开始，-90° 结束，共 181 个坐标。
  • lon_0 维度是从 0° 开始，359° 结束，共 360 个坐标。
  • 还有多组垂直维度的信息，可以参照前面查到的垂直维度信息的网页，如变量名中带 L102 的，垂直维度使用的 lv_AMSL1，其中有 3 个值，代表高于“平均海平面高度”的特定三个值。这里大家先不要头晕，只要掌握好常用的几个就好，例如 L100 代表气压层，L103 中有 2m 气温等。

另外，也可以使用以下代码，来展示变量名、变量详细信息、单位、维度等 🔰

In [4]:

# 显示文件内各个变量名及单位
for v in ds:
    print(
        "{}, {}, {}, {}".format(
            v,
            ds[v].attrs["long_name"],
            ds[v].attrs["units"],
            ds[v].shape))
TMP_P0_L1_GLL0, Temperature, K, (181, 360)
TMP_P0_L6_GLL0, Temperature, K, (181, 360)
TMP_P0_L7_GLL0, Temperature, K, (181, 360)
TMP_P0_L100_GLL0, Temperature, K, (41, 181, 360)
TMP_P0_L102_GLL0, Temperature, K, (3, 181, 360)
TMP_P0_L103_GLL0, Temperature, K, (3, 181, 360)
TMP_P0_L104_GLL0, Temperature, K, (181, 360)
TMP_P0_2L108_GLL0, Temperature, K, (181, 360)
TMP_P0_L109_GLL0, Temperature, K, (2, 181, 360)
POT_P0_L104_GLL0, Potential temperature, K, (181, 360)
DPT_P0_L103_GLL0, Dew point temperature, K, (181, 360)
APTMP_P0_L103_GLL0, Apparent temperature, K, (181, 360)
SPFH_P0_L100_GLL0, Specific humidity, kg kg-1, (41, 181, 360)
SPFH_P0_L103_GLL0, Specific humidity, kg kg-1, (2, 181, 360)
SPFH_P0_2L108_GLL0, Specific humidity, kg kg-1, (181, 360)
RH_P0_L4_GLL0, Relative humidity, %, (181, 360)
RH_P0_L100_GLL0, Relative humidity, %, (41, 181, 360)
RH_P0_L103_GLL0, Relative humidity, %, (181, 360)
RH_P0_2L104_GLL0, Relative humidity, %, (4, 181, 360)
RH_P0_L104_GLL0, Relative humidity, %, (181, 360)
RH_P0_2L108_GLL0, Relative humidity, %, (181, 360)
RH_P0_L200_GLL0, Relative humidity, %, (181, 360)
RH_P0_L204_GLL0, Relative humidity, %, (181, 360)
PWAT_P0_L200_GLL0, Precipitable water, kg m-2, (181, 360)
PRATE_P0_L1_GLL0, Precipitation rate, kg m-2 s-1, (181, 360)
SNOD_P0_L1_GLL0, Snow depth, m, (181, 360)
WEASD_P0_L1_GLL0, Water equivalent of accumulated snow depth, kg m-2, (181, 360)
CLWMR_P0_L100_GLL0, Cloud mixing ratio, kg kg-1, (22, 181, 360)
CLWMR_P0_L105_GLL0, Cloud mixing ratio, kg kg-1, (181, 360)
ICMR_P0_L100_GLL0, Ice water mixing ratio, kg kg-1, (22, 181, 360)
ICMR_P0_L105_GLL0, Ice water mixing ratio, kg kg-1, (181, 360)
RWMR_P0_L100_GLL0, Rain mixing ratio, kg kg-1, (22, 181, 360)
RWMR_P0_L105_GLL0, Rain mixing ratio, kg kg-1, (181, 360)
SNMR_P0_L100_GLL0, Snow mixing ratio, kg kg-1, (22, 181, 360)
SNMR_P0_L105_GLL0, Snow mixing ratio, kg kg-1, (181, 360)
GRLE_P0_L100_GLL0, Graupel, kg kg-1, (22, 181, 360)
GRLE_P0_L105_GLL0, Graupel, kg kg-1, (181, 360)
CPOFP_P0_L1_GLL0, Percent frozen precipitation, %, (181, 360)
CRAIN_P0_L1_GLL0, Categorical rain (yes=1; no=0), non-dim, (181, 360)
CFRZR_P0_L1_GLL0, Categorical freezing rain (yes=1; no=0), non-dim, (181, 360)
CICEP_P0_L1_GLL0, Categorical ice pellets (yes=1; no=0), non-dim, (181, 360)
CSNOW_P0_L1_GLL0, Categorical snow (yes=1; no=0), non-dim, (181, 360)
UGRD_P0_L6_GLL0, U-component of wind, m s-1, (181, 360)
UGRD_P0_L7_GLL0, U-component of wind, m s-1, (181, 360)
UGRD_P0_L100_GLL0, U-component of wind, m s-1, (41, 181, 360)
UGRD_P0_L102_GLL0, U-component of wind, m s-1, (3, 181, 360)
UGRD_P0_L103_GLL0, U-component of wind, m s-1, (7, 181, 360)
UGRD_P0_L104_GLL0, U-component of wind, m s-1, (181, 360)
UGRD_P0_2L108_GLL0, U-component of wind, m s-1, (181, 360)
UGRD_P0_L109_GLL0, U-component of wind, m s-1, (2, 181, 360)
UGRD_P0_L220_GLL0, U-component of wind, m s-1, (181, 360)
VGRD_P0_L6_GLL0, V-component of wind, m s-1, (181, 360)
VGRD_P0_L7_GLL0, V-component of wind, m s-1, (181, 360)
VGRD_P0_L100_GLL0, V-component of wind, m s-1, (41, 181, 360)
VGRD_P0_L102_GLL0, V-component of wind, m s-1, (3, 181, 360)
VGRD_P0_L103_GLL0, V-component of wind, m s-1, (7, 181, 360)
VGRD_P0_L104_GLL0, V-component of wind, m s-1, (181, 360)
VGRD_P0_2L108_GLL0, V-component of wind, m s-1, (181, 360)
VGRD_P0_L109_GLL0, V-component of wind, m s-1, (2, 181, 360)
VGRD_P0_L220_GLL0, V-component of wind, m s-1, (181, 360)
VVEL_P0_L100_GLL0, Vertical velocity (pressure), Pa s-1, (41, 181, 360)
VVEL_P0_L104_GLL0, Vertical velocity (pressure), Pa s-1, (181, 360)
DZDT_P0_L100_GLL0, Vertical velocity (geometric), m s-1, (41, 181, 360)
ABSV_P0_L100_GLL0, Absolute vorticity, s-1, (41, 181, 360)
GUST_P0_L1_GLL0, Wind speed (gust), m s-1, (181, 360)
VWSH_P0_L7_GLL0, Vertical speed shear, s-1, (181, 360)
VWSH_P0_L109_GLL0, Vertical speed shear, s-1, (2, 181, 360)
USTM_P0_2L103_GLL0, U-component storm motion, m s-1, (181, 360)
VSTM_P0_2L103_GLL0, V-component storm motion, m s-1, (181, 360)
FRICV_P0_L1_GLL0, Frictional velocity, m s-1, (181, 360)
VRATE_P0_L220_GLL0, Ventilation rate, m2 s-1, (181, 360)
PRES_P0_L1_GLL0, Pressure, Pa, (181, 360)
PRES_P0_L6_GLL0, Pressure, Pa, (181, 360)
PRES_P0_L7_GLL0, Pressure, Pa, (181, 360)
PRES_P0_L103_GLL0, Pressure, Pa, (181, 360)
PRES_P0_L109_GLL0, Pressure, Pa, (2, 181, 360)
PRMSL_P0_L101_GLL0, Pressure reduced to MSL, Pa, (181, 360)
ICAHT_P0_L6_GLL0, ICAO standard atmosphere reference height, m, (181, 360)
ICAHT_P0_L7_GLL0, ICAO standard atmosphere reference height, m, (181, 360)
HGT_P0_L1_GLL0, Geopotential height, gpm, (181, 360)
HGT_P0_L4_GLL0, Geopotential height, gpm, (181, 360)
HGT_P0_L6_GLL0, Geopotential height, gpm, (181, 360)
HGT_P0_L7_GLL0, Geopotential height, gpm, (181, 360)
HGT_P0_L100_GLL0, Geopotential height, gpm, (41, 181, 360)
HGT_P0_L109_GLL0, Geopotential height, gpm, (2, 181, 360)
HGT_P0_L204_GLL0, Geopotential height, gpm, (181, 360)
HGT_P0_L215_GLL0, Geopotential height, gpm, (181, 360)
MSLET_P0_L101_GLL0, MSLP (ETA model reduction), Pa, (181, 360)
HPBL_P0_L1_GLL0, Planetary boundary layer height, m, (181, 360)
PLPL_P0_2L108_GLL0, Pressure of level from which parcel was lifted, Pa, (181, 360)
TCDC_P0_L10_GLL0, Total cloud cover, %, (181, 360)
TCDC_P0_L100_GLL0, Total cloud cover, %, (22, 181, 360)
LCDC_P0_L214_GLL0, Low cloud cover, %, (181, 360)
MCDC_P0_L224_GLL0, Medium cloud cover, %, (181, 360)
HCDC_P0_L234_GLL0, High cloud cover, %, (181, 360)
CWAT_P0_L200_GLL0, Cloud water, kg m-2, (181, 360)
SUNSD_P0_L1_GLL0, Sunshine duration, s, (181, 360)
CAPE_P0_L1_GLL0, Convective available potential energy, J kg-1, (181, 360)
CAPE_P0_2L108_GLL0, Convective available potential energy, J kg-1, (3, 181, 360)
CIN_P0_L1_GLL0, Convective inhibition, J kg-1, (181, 360)
CIN_P0_2L108_GLL0, Convective inhibition, J kg-1, (3, 181, 360)
HLCY_P0_2L103_GLL0, Storm relative helicity, m2 s-2, (181, 360)
LFTX_P0_L1_GLL0, Surface lifted index, K, (181, 360)
4LFTX_P0_L1_GLL0, Best (4 layer) lifted index, K, (181, 360)
TOZNE_P0_L200_GLL0, Total ozone, Dobson, (181, 360)
O3MR_P0_L100_GLL0, Ozone mixing ratio, kg kg-1, (41, 181, 360)
REFD_P0_L103_GLL0, Reflectivity, dB, (2, 181, 360)
REFD_P0_L105_GLL0, Reflectivity, dB, (2, 181, 360)
REFC_P0_L10_GLL0, Composite reflectivity, dB, (181, 360)
VIS_P0_L1_GLL0, Visibility, m, (181, 360)
LAND_P0_L1_GLL0, Land cover (1=land, 0=sea), Proportion, (181, 360)
SFCR_P0_L1_GLL0, Surface roughness, m, (181, 360)
TSOIL_P0_2L106_GLL0, Soil temperature, K, (4, 181, 360)
VEG_P0_L1_GLL0, Vegetation, %, (181, 360)
SOILW_P0_2L106_GLL0, Volumetric soil moisture content, Fraction, (4, 181, 360)
CNWAT_P0_L1_GLL0, Plant canopy surface water, kg m-2, (181, 360)
WILT_P0_L1_GLL0, Wilting point, Fraction, (181, 360)
SOTYP_P0_L1_GLL0, Soil type, Code Table 4.213, (181, 360)
SOILL_P0_2L106_GLL0, Liquid volumetric soil moisture (non frozen), Proportion, (4, 181, 360)
FLDCP_P0_L1_GLL0, Field capacity, fraction, (181, 360)
HINDEX_P0_L1_GLL0, Haines index, Numeric, (181, 360)
ICEC_P0_L1_GLL0, Ice cover, Proportion, (181, 360)
ICETK_P0_L1_GLL0, Ice thickness, m, (181, 360)
ICEG_P0_L102_GLL0, Ice growth rate, m s-1, (181, 360)
ICETMP_P0_L1_GLL0, Ice temperature, K, (181, 360)
lv_DBLL11_l1, Depth below land surface, m, (4,)
lv_DBLL11_l0, Depth below land surface, m, (4,)
lv_SPDL8_l1, Level at specified pressure difference from ground to level, Pa, (3,)
lv_SPDL8_l0, Level at specified pressure difference from ground to level, Pa, (3,)
lv_SIGL5_l1, Sigma level (sigma value), none, (4,)
lv_SIGL5_l0, Sigma level (sigma value), none, (4,)

👆 注意显示的内容格式为：变量名，变量含义，单位，变量维度。

可以看到，变量名由下划线分为四个部分，第一部分是变量含义缩写，如 TMP 指的是温度，第二部分和第四部分所有变量都一样，可以忽略，第三部分是高度，如 L1 指的是地表，L6 指的是最大风速层，L7 指的是对流层顶（更多不同高度的含义见：https://rda.ucar.edu/datasets/ds083.2/#metadata/grib2_levels.html?_do=y ）。

例如：

  • TMP_P0_L100_GLL0 指的是温度（TMP）在等压面（L100）坐标的值
  • VGRD_P0_L220_GLL0 指的是经向风（VGRD）在行星边界层（L220）处的值

4. 筛选所需变量名

FNL 变量虽然比较多，但是一般不用都记住，了解你常用的几种就可以。

比如，假设我们想要读取地面或水面温度1，可以确认是变量 TMP_P0_L1_GLL0：

  1. 这里需要注意，首先要区分好到底是“地表温度”还是“地表气温”。地表温度在模型及再分析资料中，是模型中地表这一层的数据，不同的模型可能还会有差别。观测中的“地表温度”是有一套严格的标准，并且它受多种因素影响。而“地表气温”模型中多是指 2m 气温，观测中多指的是百叶箱中气温表记录的温度。↩

In [5]:

#读取TMP_P0_L1_GLL0这个变量到tt中
tt = ds.get("TMP_P0_L1_GLL0")

5. 索引你需要的数据范围

数据索引，就是我们在读取数据时，只读取我们想要的那部分。有时候变量本身比较大，而我们不需要处理全部的数据，那么在最开始读取出我们关注的那一小部分数据，可以提高代码运行速度。本文使用的 xarray 提供了三种方式：下标索引、坐标值索引和标签索引。由于标签索引使用的比较少，以下主要介绍下标索引和坐标值索引。

下标索引指的是针对数据本身的结构，如

In [6]:

tt1 = ds.get("TMP_P0_L1_GLL0")[0:4,11:15]
tt1

In [7]:

tt2 = ds.get("TMP_P0_L1_GLL0").loc[35:31,100:115]
tt2

上述代码中：

  • tt1 指的是读取地表温度的纬度中第 0 到 4、经度中第 11 到 15 的范围。通过数据概况可知，这代表纬度为北纬 87°-90°、东经 11°-14° 的范围。
  • tt2 读取的则是纬度为北纬 31°-35°、东经 100°-115° 的范围。注意这里的顺序要写为 35:31，📌因为源数据的纬度是从北到南的，所以索引数据也需要按照从北到南来写。

再来一个读取 2022 年 1 月 1 日 0 时北半球寒带地球地表以上 2m 气温的范例：

In [8]:

ds_camp = xr.open_dataset("/home/mw/input/ncep4389/fnl_20220101_06_00.grib2", engine="pynio") # 读取2022010100FNL数据
ds_camp

在 FNL 数据对垂直层的介绍中可知，L103 指的是“Specified height level above ground”。在上方输出中，点击 lv_HTGL2 的说明（show/hide data repr），可以看到有 array([ 2., 80., 100.], dtype=float32) 的式样，说明这三个维度分别代表 2m、80m、100m 的高度。

因此：

  • 2 米温度对应要选择 TMP_P0_L103_GLL0 的变量中高度第 0 维的部分，所以你可以通过代码 [0,....] 来索引。
  • 北半球寒带地区，维度从北极到北纬 66°34′，所以你可以通过代码 .loc[:66+34/60:] 来索引。

为了向你演示过程，方便你理解地表以上 2m 温度数据，这里我们分两次索引

（实际你可以直接使用 ds_camp["TMP_P0_L103_GLL0"][0,...].loc[:66+34/66,:] 代码直接完成索引）。

In [9]:

t_camp = ds_camp["TMP_P0_L103_GLL0"][0,...] #索引地表以上2m温度
t_camp #查看地表以上2m温度

In [10]:

t_north = t_camp.loc[:66+34/60,:] #索引北半球寒带数据
t_north
FEATURED TAGS
科研 博客 漫谈 satellite models Emission Inventory plots python matplotlib emissions
FRIENDS