使用 R 语言处理 netCDF 格式的数据

应 @gwy 小伙伴的定制，特推出使用 R 语言处理 netCDF 格式数据的教程，本教程包含两次可，主要是带领大家了解 netCDF 格式数据的结构、属性提取、空间插值与可视化等操作。

netCDF 的介绍

NetCDF(network Common Data Form) 网络通用数据格式是由美国大学大气研究协会（University Corporation for Atmospheric Research，UCAR)的 Unidata 项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF 广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。用户可以借助多种方式方便地管理和操作 NetCDF 数据集。（摘抄自百度百科）

本课程包含如下内容：

1. 将 netCDF 格式的数据读入 R；
2. netCDF 格式数据的解读；
3. 如何将 netCDF 数据转成数据框；
4. 地理数据可视化；
5. 计算中国各县的平均地表气压；
6. 空间插值的使用。

课程材料预览

节选自「使用R语言处理 netCDF 格式的数据（二）」

继续之前的使用 R 语言处理 netCDF 格式的数据。上次我们使用的是从 NASA 下载的温度数据演示的，但是那个数据存在很多缺失，缺失的原因我也在视频课上说了，所以这次就换了一个，也是从 https://disc.gsfc.nasa.gov/datasets/M2TMNXAER_5.12.4/summary 下载的，不过是地表气压（surface air pressure）数据。

本次课程首先回顾了上次课的内容，介绍了如何使用 R 读取 netCDF 数据以及如何从读取得到的 ncdf4 对象中提取自己需要的数据。由于 NASA 提供的数据很粗（格点不够密集），如果直接用来计算每个县的平均地表气压会有很多县的数据是缺失的（这些县的面积小，没有任何格点落在县域里），所以可以先进行空间插值的到更为细密的格点，然后再使用插值后的数据计算各个县的平均地表气压。

首先加载我们需要的两个包：

library(ncdf4)
library(tidyverse)

打开本次做展示的 netCDF 文件:

nc_open('MERRA2_400.inst6_3d_ana_Np.20200430.SUB.nc') -> ncin
print(ncin)

打印的信息很长，我就不 po 出来了，其他的信息和上次课中讲的基本一样，这次我们选择 PS （surface air pressure）数据进行处理，这个数据只有三个维度：

float PS[lon,lat,time]   (Chunking: [144,91,1])  (Compression: shuffle,level 1)
    standard_name: surface_air_pressure
    long_name: Surface pressure
    units: Pa
    _FillValue: 999999986991104
    missing_value: 999999986991104
    fmissing_value: 999999986991104
    vmax: 999999986991104
    vmin: -999999986991104

下面我们按照和上次课一样的思路把 PS 数据提取为数据框：

# 获取经度和纬度
ncvar_get(ncin, "lon") -> lon
ncvar_get(ncin, "lat") -> lat

...此处隐藏部分代码...

pssf
#> Simple feature collection with 17545 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 60 ymin: 0 xmax: 150 ymax: 60
#> CRS:            EPSG:4326
#> # A tibble: 17,545 x 2
#>       ps    geometry
#>  * <dbl> <POINT [°]>
#>  1 1010.      (60 0)
#>  2 1010.    (60 0.5)
#>  3 1010.      (60 1)
#>  4 1010.    (60 1.5)
#>  5 1010.      (60 2)
#>  6 1010.    (60 2.5)
#>  7 1010.      (60 3)
#>  8 1009.    (60 3.5)
#>  9 1010.      (60 4)
#> 10 1010.    (60 4.5)
#> # … with 17,535 more rows

现在的到的格点数量是 17545 个。

# 再读取一个中国地图数据
read_sf('provmap2019new.geojson') -> cn

绘图展示地表气压数据：

ggplot() +
  geom_sf(data = pssf,
          aes(color = ps), size = 0.01) +
  geom_sf(data = cn, size = 0.3, fill = NA) +
  labs(color = "地表气压（hPa）") +
  theme_ipsum(base_family = cnfont)

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言处理 netCDF 格式的数据