在R中从netCDF中提取时间序列
我在1998年1月1日至1998年12月31date间使用TRMM_3B42_Daily产品创build了这个文件 。 这是我在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) 在这个链接之后，我尝试提取这些值来绘制一个时间序列，但是它似乎是对两个单元格的值进行平均，而不是仅仅提取单个单元格的值。 我该如何解决？ 有没有办法创build一个循环，以便提取不同单元格的值？ （在这种情况下，它将是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文件中的date在前面有一个X，值是水平显示，而不是垂直显示。 任何帮助将不胜感激！！ 谢谢
点数据的提取可以通过创build包含要提取数据的点的SpatialPoints对象，然后执行extract操作来轻松完成。 关于其他主题：添加“X”是因为列名不能以数字开头，所以添加了一个字符。 提取一些转置后，可以很容易地改变水平sorting
例如，这应该工作（它也解决了“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`)