最近有小伙伴问到了关于如何处理从 https://disc.gsfc.nasa.gov/datasets/M2I1NXLFO_5.12.4/summary 页面下载到的近地表的气压、比湿、风速和气温的数据。直接下载到的数据为 NetCDF 格式，虽然之前也讲解过类似的栅格数据转面板课程，不过感觉之前的代码都不够简洁。今天的课程中我们来一起看看如何高效的处理这份数据。
该数据包含如下变量：
其中比湿（specific humidity）是指在一团湿空气中，水汽的质量与该团空气总质量（水汽质量加上干空气质量）的比值。如果湿空气与外界无质量交换，并且没有相变发生，那么比湿会保持不变。比湿通常以g/g或g/kg为单位来表示，通常大气中的比湿都小于40g/kg，是记录大气水汽状况的一个重要指标。
下面我们一起来看下这个数据该如何下载和处理。
从 NASA 下载 MERRA-2 数据
首先你需要注册一个 NASA 的账号，这里就不再演示了：https://disc.gsfc.nasa.gov/datasets/M2I1NXLFO_5.12.4/summary
点击右下角的 Subset/Get Data 就可以选择自己需要的数据下载了：
- Download Method 选择 Get File Subsets using OPeNDAP；
- Refine Date Range 选择 2023-01-01～2023-12-31；
- Refine Region 选择 73.502,3.837,135.096,53.564，这个是中国的经纬度范围：
library(sf)
sf::read_sf("2021行政区划/省.shp") -> prov
st_bbox(prov) %>%
round(3) %>%
paste0(collapse = ",")
|
#> [1] "73.502,3.837,135.096,53.564"
|
然后变量都选上：
最后点击 Get Data 就可以获取数据下载链接了：
点击 Download links list 下载所有的数据链接，就是附件中的 subset_M2I1NXLFO_5.12.4_20240515_082340_.txt 文件，不过需要注意里面的数据链接有效期仅为两天（大家需要自己重新生成这个链接）。
然后我们就可以使用 NASA 上面介绍的方法下载这些数据了：
由于这些数据需要登录验证才能下载，所以大家可以一个个把链接复制到浏览器的地址栏回车下载，也可以使用官网介绍的一些方法批量下载：https://disc.gsfc.nasa.gov/data-access
这里我演示下 Mac 上使用 curl 批量下载的方法：
这里演示的方法仅适用于 Mac。由于我手上没有 Windows 电脑，所以如果大家使用 Windows 话可以对着上述网页的介绍自己尝试下。
打开终端（terminal）依次运行下面的代码（记得先把 替换成自己的 NASA 账户名， 替换成自己的账户密码）：
cd ~
touch .netrc
echo "machine urs.earthdata.nasa.gov login <uid> password <password>" >> .netrc
chmod 0600 .netrc
touch .urs_cookies
|
然后我们在工作目录下面新建一个 data 文件夹，把前面下载到的 subset_M2I1NXLFO_5.12.4_20240515_082340_.txt 文件放进去，然后在终端运行下面的代码（需要把 “data 文件的完整路径” 替换成 data 文件的完整路径）：
cd "/Users/ac/Desktop/使用 R 语言处理 Merra2 数据获取各省市区县的比湿、降水量、风速和气压数据/data/"
cat subset_M2I1NXLFO_5.12.4_20240515_082340_.txt | tr -d '\r' | xargs -n 1 curl --globoff -LJO -n -c ~/.urs_cookies -b ~/.urs_cookies
|
注意这里我加了 --globoff 参数，要不然会提示说链接无效。
然后就开始下载数据了。
一个个的下载会很慢，所以也可以把下载链接拆分成多个文件，然后打开多个终端窗口同时下载：
library(tidyverse)
read_lines("data/subset_M2I1NXLFO_5.12.4_20240515_082340_.txt") %>%
as_tibble() %>%
rename(X1 = value) -> df
df
df %>%
slice(-1) -> df
df %>%
mutate(filename = str_match(X1, "\\d{8}")[,1]) %>%
anti_join(
fs::dir_ls("data") %>%
as.character() %>%
as_tibble() %>%
mutate(filename = str_match(value, "\\d{8}")[,1])
) %>%
select(X1) -> df
lapply(seq(1, nrow(df), by = 75), function(x){
df %>%
slice(x:(x+74)) %>%
write_csv(paste0("data/downlist", x, ".txt"), col_names = F, quote = "none")
}) -> res
for (x in seq(1, nrow(df), by = 75)) {
message("cd '/Users/ac/Desktop/使用 R 语言处理 Merra2 数据获取各省市区县的比湿、降水量、风速和气压数据/data/'")
paste0("cat downlist", x, ".txt | tr -d '\\r' | xargs -n 1 curl --globoff -LJO -n -c ~/.urs_cookies -b ~/.urs_cookies") %>% message()
message(" ")
}
|
由于 curl 下载经常遇到各种网络原因中断，所以需要反复下载、剔除下载错误的。
这样我们就下载到了 2023 年的数据。由于时间的关系，我只下载了 10 天的。
读取 nc 文件
使用 terra 包就可以读取 NetCDF 文件了，不过经常无法识别坐标系，所以需要根据实际情况来设定，可能会需要进行旋转(trans)、镜像变换(flip)、设置 crs 之类的。
library(tidyverse)
library(sf)
library(terra)
|
注意我这里使用的是 1.7.71 版本的 terra 包。
# 读取 nc 文件
rast("data/MERRA2_400.inst1_2d_lfo_Nx.20230219.nc4.nc4") -> rst
plot(rst$TLML_24)
|
如果发现读取的结果不对，可以考虑使用 trans() 和 flip() 进行旋转和镜像变换。
如果发现没有识别到 extent，可以用下面的代码设置范围：
ext(rst) <- c(73.4375, 135.3125, 3.75, 53.75)
|
如果发现没有识别到 crs，可以用下面的代码设置坐标系：
crs(rst) <- "+proj=longlat +datum=WGS84 +no_defs"
|
由于这个数据是个小时数据，所以可以通过平均得到日度的，一个直观的方法是加总再平均：
(rst$HLML_1 + rst$HLML_2 + rst$HLML_3 +
rst$HLML_4 + rst$HLML_5 + rst$HLML_6 +
rst$HLML_7 + rst$HLML_8 + rst$HLML_9 +
rst$HLML_10 + rst$HLML_11 + rst$HLML_12 +
rst$HLML_13 + rst$HLML_14 + rst$HLML_15 +
rst$HLML_16 + rst$HLML_17 + rst$HLML_18 +
rst$HLML_19 + rst$HLML_20 + rst$HLML_21 +
rst$HLML_22 + rst$HLML_23 + rst$HLML_24)/24 -> newrst
plot(newrst)
|
由于 HLML 指标是 rst 的 1～24 层，所以也可以用下面的代码进行平均：
app(rst[[1:24]], fun = mean) -> newrst
plot(newrst)
|
由于这个文件含有多个指标，可以一并处理：
rast(
list(
app(rst[[1:24]], fun = mean),
app(rst[[25:48]], fun = mean),
app(rst[[49:72]], fun = mean),
app(rst[[73:96]], fun = mean),
app(rst[[97:120]], fun = mean)
)
) %>%
`names<-`(varnames(rst)) -> newrst
newrst
|
#> class : SpatRaster
#> dimensions : 100, 99, 5 (nrow, ncol, nlyr)
#> resolution : 0.625, 0.5 (x, y)
#> extent : 73.4375, 135.3125, 3.75, 53.75 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> source(s) : memory
#> names : HLML, PS, QLML, SPEEDLML, TLML
#> min values : 55.44875, 51040.24, 0.0004833084, 1.041553, 250.6861
#> max values : 67.51266, 102821.22, 0.0198902683, 18.710755, 303.1523
|
为了方便使用，我们可以把这个对象再保存成 tif 文件：
newrst %>%
writeRaster("temp.tif", overwrite = T)
|
栅格数据转面板
我们经常会需要把栅格数据汇总成各省市区县的面板来使用。
读取省市区县的 shp 文件：
read_sf("2021行政区划/省.shp") %>%
select(-contains("类型")) %>%
st_transform(4326) -> prov
read_sf("2021行政区划/市.shp") %>%
select(-contains("类型")) %>%
st_transform(4326) -> city
read_sf("2021行政区划/县.shp") %>%
select(-contains("类型")) %>%
st_transform(4326) -> county
|
使用 terra::extract() 提取：
dir.create("restemp")
y <- str_extract("data/MERRA2_400.inst1_2d_lfo_Nx.20230219.nc4.nc4", "\\d{8}")
terra::extract(newrst, vect(prov), fun = mean, na.rm = T) %>%
as_tibble() -> dfres1
terra::extract(newrst, vect(city), fun = mean, na.rm = T) %>%
as_tibble() -> dfres2
terra::extract(newrst, vect(county), fun = mean, na.rm = T) %>%
as_tibble() -> dfres3
prov %>%
st_drop_geometry() %>%
bind_cols(dfres1) %>%
mutate(date = y) %>%
write_rds(paste0("restemp/", y, "_prov.rds"))
city %>%
st_drop_geometry() %>%
bind_cols(dfres2) %>%
mutate(date = y) %>%
write_rds(paste0("restemp/", y, "_city.rds"))
county %>%
st_drop_geometry() %>%
bind_cols(dfres3) %>%
mutate(date = y) %>%
write_rds(paste0("restemp/", y, "_county.rds"))
|
循环处理所有的
如果文件较少、每个文件只有一个指标，可以使用 rast() 直接读取所有的文件一起处理，代码会简洁很多，不过因为这里每个文件都有 5 个指标，所以混在一起会难以区分。
fs::dir_ls("data", regexp = "nc4$") -> ls
rast(ls) -> rst
rst
|
#> class : SpatRaster
#> dimensions : 100, 99, 1200 (nrow, ncol, nlyr)
#> resolution : 0.625, 0.5 (x, y)
#> extent : 73.4375, 135.3125, 3.75, 53.75 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84
#> sources : MERRA2_400.inst1_2d_lfo_Nx.20230210.nc4.nc4:HLML (24 layers)
#> MERRA2_400.inst1_2d_lfo_Nx.20230210.nc4.nc4:PS (24 layers)
#> MERRA2_400.inst1_2d_lfo_Nx.20230210.nc4.nc4:QLML (24 layers)
#> ... and 47 more source(s)
#> varnames : HLML (surface_layer_height)
#> PS (surface_pressure)
#> QLML (surface_specific_humidity)
#> ...
#> names : HLML_1, HLML_2, HLML_3, HLML_4, HLML_5, HLML_6, ...
#> unit : m, m, m, m, m, m, ...
#> time : 2023-02-10 00:00:00 to 2023-02-19 23:00:00 UTC
|
因此我们这里还是不要这样处理了。
循环文件更好：
dir.create("tif")
dir.create("res")
lapply(ls, function(x){
print(x)
rast(x) -> rst
rast(
list(
app(rst[[1:24]], fun = mean),
app(rst[[25:48]], fun = mean),
app(rst[[49:72]], fun = mean),
app(rst[[73:96]], fun = mean),
app(rst[[97:120]], fun = mean)
)
) %>%
`names<-`(varnames(rst)) -> newrst
str_extract(x, "\\d{8}") -> y
newrst %>%
writeRaster(paste0("tif/", y, ".tif"))
terra::extract(newrst, vect(prov), fun = mean, na.rm = T) %>%
as_tibble() %>%
mutate(date = y) %>%
write_rds(paste0("res/", y, "_prov.rds"))
terra::extract(newrst, vect(city), fun = mean, na.rm = T) %>%
as_tibble() %>%
mutate(date = y) %>%
write_rds(paste0("res/", y, "_city.rds"))
terra::extract(newrst, vect(county), fun = mean, na.rm = T) %>%
as_tibble() %>%
mutate(date = y) %>%
write_rds(paste0("res/", y, "_county.rds")) -> dfres3
}) -> res
|
这里如果文件很多，可以用多线程：
library(parallel)
detectCores()
|
makeCluster(8) -> cl
clusterEvalQ(cl, ({
library(tidyverse)
library(terra)
library(sf)
}))
clusterExport(cl, "prov")
clusterExport(cl, "city")
clusterExport(cl, "county")
|
parLapply(cl, ls, function(x){
str_extract(x, "\\d{8}") -> y
if(!file.exists(paste0("res/", y, "_county.rds"))) {
rast(x) -> rst
rast(
list(
app(rst[[1:24]], fun = mean),
app(rst[[25:48]], fun = mean),
app(rst[[49:72]], fun = mean),
app(rst[[73:96]], fun = mean),
app(rst[[97:120]], fun = mean)
)
) %>%
`names<-`(varnames(rst)) -> newrst
newrst %>%
writeRaster(paste0("tif/", y, ".tif"), overwrite = T)
terra::extract(newrst, vect(prov), fun = mean, na.rm = T) %>%
as_tibble() %>%
mutate(date = y) %>%
write_rds(paste0("res/", y, "_prov.rds"))
terra::extract(newrst, vect(city), fun = mean, na.rm = T) %>%
as_tibble() %>%
mutate(date = y) %>%
write_rds(paste0("res/", y, "_city.rds"))
terra::extract(newrst, vect(county), fun = mean, na.rm = T) %>%
as_tibble() %>%
mutate(date = y) %>%
write_rds(paste0("res/", y, "_county.rds"))
}
}) -> res
|
合并处理结果
最后就可以分别合并省市区县的结果了：
fs::dir_ls("res", regexp = "prov") %>%
parLapply(cl, ., read_rds) %>%
bind_rows() -> df1
prov %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
right_join(df1) %>%
select(-ID) %>%
select(1:2, date, everything()) %>%
mutate(date = ymd(date))
|
#> # A tibble: 340 × 8
#> 省 省代码 date HLML PS QLML SPEEDLML TLML
#> <chr> <dbl> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 安徽省 340000 2023-02-10 61.8 100755. 0.00455 4.54 279.
#> 2 安徽省 340000 2023-02-11 62.3 100413. 0.00558 5.07 281.
#> 3 安徽省 340000 2023-02-12 62.2 100542. 0.00565 7.53 280.
#> 4 安徽省 340000 2023-02-13 61.0 101502. 0.00340 7.45 275.
#> 5 安徽省 340000 2023-02-14 61.1 101694. 0.00243 4.48 276.
#> 6 安徽省 340000 2023-02-15 61.4 101867. 0.00260 3.99 277.
#> 7 安徽省 340000 2023-02-16 62.1 101380. 0.00350 3.70 280.
#> 8 安徽省 340000 2023-02-17 62.9 100764. 0.00494 3.89 283.
#> 9 安徽省 340000 2023-02-18 62.8 100991. 0.00548 6.07 283.
#> 10 安徽省 340000 2023-02-19 62.2 101443. 0.00310 4.56 281.
#> # ℹ 330 more rows
|
fs::dir_ls("res", regexp = "city") %>%
parLapply(cl, ., read_rds) %>%
bind_rows() -> df2
city %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
right_join(df2) %>%
select(-ID) %>%
select(1:4, date, everything()) %>%
mutate(date = ymd(date))
|
#> # A tibble: 3,710 × 10
#> 省 省代码 市 市代码 date HLML PS QLML SPEEDLML TLML
#> <chr> <dbl> <chr> <dbl> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 安徽省 340000 安庆市 340800 2023-02-10 61.8 99519. 0.00481 4.32 278.
#> 2 安徽省 340000 安庆市 340800 2023-02-11 62.5 99124. 0.00592 4.55 282.
#> 3 安徽省 340000 安庆市 340800 2023-02-12 62.5 99238. 0.00614 6.62 281.
#> 4 安徽省 340000 安庆市 340800 2023-02-13 61.1 100210. 0.00394 8.21 275.
#> 5 安徽省 340000 安庆市 340800 2023-02-14 61.1 100450. 0.00274 5.28 276.
#> 6 安徽省 340000 安庆市 340800 2023-02-15 61.5 100644. 0.00280 4.11 278.
#> 7 安徽省 340000 安庆市 340800 2023-02-16 62.2 100192. 0.00343 2.29 281.
#> 8 安徽省 340000 安庆市 340800 2023-02-17 63.0 99623. 0.00532 4.28 284.
#> 9 安徽省 340000 安庆市 340800 2023-02-18 63.0 99821. 0.00624 5.54 284.
#> 10 安徽省 340000 安庆市 340800 2023-02-19 62.4 100242. 0.00352 5.96 281.
#> # ℹ 3,700 more rows
|
fs::dir_ls("res", regexp = "county") %>%
parLapply(cl, ., read_rds) %>%
bind_rows() -> df3
county %>%
st_drop_geometry() %>%
mutate(ID = row_number()) %>%
right_join(df3) %>%
select(-ID) %>%
select(1:6, date, everything()) %>%
mutate(date = ymd(date))
|
#> # A tibble: 28,770 × 12
#> 省 省代码 市 市代码 县 县代码 date HLML PS QLML
#> <chr> <dbl> <chr> <dbl> <chr> <dbl> <date> <dbl> <dbl> <dbl>
#> 1 安徽省 340000 安庆市 340800 大观区 340803 2023-02-10 62.0 100939. 0.00501
#> 2 安徽省 340000 安庆市 340800 大观区 340803 2023-02-11 62.6 100529. 0.00598
#> 3 安徽省 340000 安庆市 340800 大观区 340803 2023-02-12 62.6 100631. 0.00630
#> 4 安徽省 340000 安庆市 340800 大观区 340803 2023-02-13 61.2 101634. 0.00414
#> 5 安徽省 340000 安庆市 340800 大观区 340803 2023-02-14 61.2 101890. 0.00280
#> 6 安徽省 340000 安庆市 340800 大观区 340803 2023-02-15 61.6 102087. 0.00290
#> 7 安徽省 340000 安庆市 340800 大观区 340803 2023-02-16 62.4 101619. 0.00335
#> 8 安徽省 340000 安庆市 340800 大观区 340803 2023-02-17 63.2 101038. 0.00530
#> 9 安徽省 340000 安庆市 340800 大观区 340803 2023-02-18 63.2 101207. 0.00639
#> 10 安徽省 340000 安庆市 340800 大观区 340803 2023-02-19 62.5 101652. 0.00381
#> # ℹ 28,760 more rows
#> # ℹ 2 more variables: SPEEDLML <dbl>, TLML <dbl>
|
这样就处理好了～
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言处理 Merra2 数据获取各省市区县的比湿、降水量、风速和气压数据
评论