2017-02-07 17 views
1

私はこれを作成しましたfile TRMM_3B42_Daily製品を使用して1998-01-01〜1998-12-31を作成しました。これは私がRで使用されるスクリプトです:netCDFから時系列を抽出する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後、私は時系列をプロットするために、値を抽出しようとしたが、私が2つのセルの値を平均化するだけではなく、Aの値を抽出していそうです単細胞。これをどうやって解決するのですか?ループを作成して、異なるセルの値を抽出する方法はありますか? (この場合には、13×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を有し、値は水平の代わりに垂直に示されていることをことをどこ。どんな助けでも大歓迎です!!ありがとう

答えて

1

ポイントデータの抽出は、extract操作のあとに、データを抽出するポイントを含む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`)