从 R 中的 netCDF 中提取时间序列

Extract time series from netCDF in R

我创造了这个file 在 1998-01-01 到 1998-12-31 期间使用 TRMM_3B42_Daily 产品。这是我在 R:

中使用的脚本
lon=seq(-91.875,-86.875,by= 0.25)
lat=seq(13.875,16.875,by= 0.25)

x_dim <- ncdim_def( "lon", "degrees_east", lon, create_dimvar=TRUE)
y_dim <- ncdim_def( "lat", "degrees_north", lat, create_dimvar=TRUE)
t_dim <- ncdim_def( "time", "days since 1997-12-31 12:00:00.0 -0:00", 1:365, unlim=FALSE)
mv=9999.900390625 
precipitation_var <- ncvar_def("precipitation", "mm", list(y_dim,x_dim,t_dim), mv)


nrow = 13 
ncol = 21 

NA.matrix=matrix(rep(NA,nrow*ncol)) 

precip=array(NA.matrix,c(nrow,ncol, 1))
for (i in 1:length(test01)){precip_nc=nc_open(test01[i])
precip_get_nc=ncvar_get(precip_nc,"precipitation") 
precip=abind(precip,precip_get_nc)}

precip=precip[,,-1]  

PRECIPITATION_nc = nc_create("PRECIPITATION_1998.nc", precipitation_var)

precipitation_nc_put=ncvar_put (PRECIPITATION_nc, precipitation_var, precip)

nc_close(PRECIPITATION_nc)

在此之后 link 我尝试提取值以绘制时间序列,但似乎我是在对两个单元格的值进行平均,而不是仅仅提取单个单元格的值。我该如何解决?有没有办法创建一个循环,以便它提取不同单元格的值? (在本例中为 13 x 21 = 273)

b <- brick('PRECIPITATION_1998.nc')
be <- crop(b, extent(13.875, 14.125, -91.875,-91.625))
a <- aggregate(be, dim(be)[2:1], na.rm=TRUE)
v <- values(a)
write.csv(v, 'precip.csv', row.names=FALSE)

我还发现了另外两个问题,即 excel 文件中的日期前面有一个 X,并且值是水平显示的而不是垂直显示的。任何帮助将不胜感激!!谢谢

通过创建包含要从中提取数据的点的 SpatialPoints 对象,然后执行 extract 操作,可以轻松完成点数据的提取。 关于其他主题: 添加"X"是因为列名不能以数字开头,所以添加了一个字符。提取后可以通过一些转置轻松更改水平顺序

例如,这应该有效(它还解决了 "X" 的问题并将格式更改为 "column like"):

library(raster)
library(stringr)
library(lubridate)
library(tidyverse)

b <- brick('/home/lb/Temp/buttami/PRECIPITATION_1998.nc')
lon = c(-91.875,-91.625)  # Array of x coordinates
lat <- c(13.875, 14.125)  # Array of y coordinates
points <- SpatialPoints(cbind(lat,lon)), # Build a spPoints object

# Etract and tidy
points_data <- b %>% 
  raster::extract(points, df = T) %>% 
  gather(date, value, -ID) %>% 
  spread(ID, value) %>%   # Can be skipped if you want a "long" table
  mutate(date = ymd(str_sub(names(b),2))) %>% 
  as_tibble()

points_data 

# A tibble: 365 × 3
         date   `1`   `2`
       <date> <dbl> <dbl>
1  1998-01-01     0     0
2  1998-01-02     0     0
3  1998-01-03     0     0
4  1998-01-04     0     0
5  1998-01-05     0     0
6  1998-01-06     0     0
7  1998-01-07     0     0
8  1998-01-08     0     0
9  1998-01-09     0     0
10 1998-01-10     0     0
# ... with 355 more rows

plot(points_data$date,points_data$`1`)