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

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

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

视频讲解

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

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

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

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

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

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

1
2
3
4
5
6
7
8
9
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 数据提取为数据框:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

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

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

绘图展示地表气压数据:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

可以看到格点的间距很大,所以为了计算各个县的平均,我们需要进行插值的到更细密的网格。

空间插值可以使用 gstat,这个包提供了多种空间插值的方法,今天我们使用 IDW 法(反向距离加权),不过在此之前我们可以先从刚刚的到的 pssf 数据中提取位于中国内部的格点(这样可以减少插值时的计算量):

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#> Simple feature collection with 3088 features and 2 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 73.75 ymin: 18.5 xmax: 134.375 ymax: 53.5
#> CRS: EPSG:4326
#> # A tibble: 3,088 x 3
#> ps 省份 geometry
#> <dbl> <chr> <POINT [°]>
#> 1 887. 广西壮族自治区 (105 24.5)
#> 2 913. 广西壮族自治区 (105.625 24.5)
#> 3 930. 广西壮族自治区 (106.25 23)
#> 4 924. 广西壮族自治区 (106.25 23.5)
#> 5 950. 广西壮族自治区 (106.25 24)
#> 6 917. 广西壮族自治区 (106.25 24.5)
#> 7 980. 广西壮族自治区 (106.875 22)
#> 8 981. 广西壮族自治区 (106.875 22.5)
#> 9 954. 广西壮族自治区 (106.875 23)
#> 10 957. 广西壮族自治区 (106.875 23.5)
#> # … with 3,078 more rows

注意到这个时候 pssf 只剩下 3088 个观测值了,而中国有 2800 多个县级行政单位,显然会有很多县没有分配到任何观测值。

下面使用 gstat 包的 idw 函数应用 IDW 插值。这个函数实际上可以传入 sf 对象,但是我测试了下,使用 sf 对象进行插值会很卡很慢,所以我还是先把 sf 对象转换为 sp 对象进行运算:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

打印出来看一下我们的结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
psidwsf
#> Simple feature collection with 6268355 features and 1 field
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 73.49008 ymin: 3.40848 xmax: 135.0863 ymax: 53.55999
#> CRS: EPSG:4326
#> # A tibble: 6,268,355 x 2
#> ps geometry
#> <dbl> <POINT [°]>
#> 1 854. (111.7721 3.40848)
#> 2 854. (111.8029 3.40848)
#> 3 854. (111.8029 3.433556)
#> 4 854. (111.7721 3.433556)
#> 5 854. (111.7721 3.40848)
#> 6 854. (111.8029 3.40848)
#> 7 854. (111.8337 3.40848)
#> 8 854. (111.8337 3.433556)
#> 9 854. (111.8029 3.433556)
#> 10 854. (111.8029 3.40848)
#> # … with 6,268,345 more rows

得到了将近 627 万个格点,应该够了,接下来再读入县级地图的数据以及计算各个县的平均地表气压:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

这个时候是否还有缺失的呢?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
countyps %>% 
dplyr::filter(is.na(ps))
#> Simple feature collection with 12 features and 2 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: 116.9403 ymin: 22.61624 xmax: 123.7016 ymax: 39.3054
#> CRS: 4326
#> First 10 features:
#> ps geometry code
#> 1 NA MULTIPOLYGON (((122.7201 39... 210224
#> 2 NA MULTIPOLYGON (((122.0684 30... 330922
#> 3 NA MULTIPOLYGON (((118.1652 24... 350527
#> 4 NA MULTIPOLYGON (((120.6868 37... 370634
#> 5 NA MULTIPOLYGON (((117.2862 23... 440523
#> 6 NA MULTIPOLYGON (((123.5444 25... 710024
#> 7 NA MULTIPOLYGON (((119.4435 23... 710036
#> 8 NA MULTIPOLYGON (((121.5329 25... 710109
#> 9 NA MULTIPOLYGON (((121.5223 25... 710111
#> 10 NA MULTIPOLYGON (((120.3167 22... 710206

呵!还有 12 个,随它们去吧!

最后展示下插值后的成果:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

这样我们就计算出来 4 月 30 日 0 点的时候的中国各县平均地表气压数据了。只要通过循环即可获取所有天的。我就不再赘述了。

知识星球附件链接:https://t.zsxq.com/2rjYJEM

#

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×