从 RasterBrick 中提取特定坐标的数据

Extracting data for certain coordinates from RasterBrick

来自 this file 我想提取特定坐标的所有值。

所以我创建了一个 RasterBrick:

library(raster)
library(R.utils)

rr <- gunzip("rr_0.25deg_reg_v12.0.nc.gz")
rr <- brick(rr)

> rr
class       : RasterBrick 
dimensions  : 201, 464, 93264, 23922  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : -40.5, 75.5, 25.25, 75.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\...\rr_0.25deg_reg_v12.0.nc 
names       : X1950.01.01, X1950.01.02, X1950.01.03, X1950.01.04, X1950.01.05, X1950.01.06, X1950.01.07, X1950.01.08, X1950.01.09, X1950.01.10, X1950.01.11, X1950.01.12, X1950.01.13, X1950.01.14, X1950.01.15, ... 
Date        : 1950-01-01, 2015-06-30 (min, max)
varname     : rr 

现在我想提取某些坐标的所有值(从 1950-01-01 到 2015-06-30 的每日降水量),例如经度 9.258157 纬度 47.70842

我试过 extract 是这样的:

raster::extract(rr, cbind(9.258157, 47.70842), df = TRUE)

导致:

'data.frame':   1 obs. of  23923 variables:
 $ ID         : num 1
 $ X1950.01.01: num 0
 $ X1950.01.02: num 9.3
 $ X1950.01.03: num 4.2
 $ X1950.01.04: num 11.2
 $ X1950.01.05: num 0
 $ X1950.01.06: num 0
 $ X1950.01.07: num 0
 $ X1950.01.08: num 0
 $ X1950.01.09: num 0
 $ X1950.01.10: num 0
 ...

但我想要一个包含两列(日期和值)的数据框作为输出,所以我做了:

rr.plot <- t(raster::extract(rr, cbind(9.258157, 47.70842), df = TRUE))
rr.plot <- data.frame(Date = as.Date(substr(rownames(rr.plot), 2, 11), format = "%Y.%m.%d"), Prcp = rr.plot[,1])[-1,]
rownames(rr.plot) <- 1:nrow(rr.plot)

> head(rr.plot)
        Date Prcp
1 1950-01-01  0.0
2 1950-01-02  9.3
3 1950-01-03  4.2
4 1950-01-04 11.2
5 1950-01-05  0.0
6 1950-01-06  0.0

但我怀疑这是最好的方法,一定有更直接的解决方案吗?

使用tidyr::gather怎么样?

xy <- raster::extract(rr, cbind(9.258157, 47.70842), df = TRUE)[, 1:10]

library(tidyr)
xy
  ID X1950.01.01 X1950.01.02 X1950.01.03 X1950.01.04 X1950.01.05 X1950.01.06 X1950.01.07 X1950.01.08 X1950.01.09
1  1           0         9.3         4.2        11.2           0           0           0           0           0

gather(xy[, -1], key = datum, value = Prcp)
        datum Prcp
1 X1950.01.01  0.0
2 X1950.01.02  9.3
3 X1950.01.03  4.2
4 X1950.01.04 11.2
5 X1950.01.05  0.0
6 X1950.01.06  0.0
7 X1950.01.07  0.0
8 X1950.01.08  0.0
9 X1950.01.09  0.0