从R中的netCDF提取时间序列

Question

我在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)

通过这个链接，我试图提取值来绘制时间序列，但似乎我正在平均两个单元格的值，而不仅仅是提取单个单元格的值。我该如何解决这个问题？有没有办法创建一个循环，以便它提取不同单元格的值？（在这种情况下，它将是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，并且值水平显示而不是垂直显示。任何帮助将不胜感激！谢谢

Answer 1

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

例如，这应该有效（它还解决了“X”问题并将格式更改为“类似列”）：

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`)