2017-11-26 6 views
-1

私はこのベクトル(時系列)計算コードを再現しようとしている:ラスター計算コードにr - この時系列計算をラスタ計算に変換するにはどうすればよいですか?

gamma.parameters<- fitdistr(may_baseline_3months[may_baseline_3months>0],"gamma")

このコードは元々、ベクトル(時系列)may_baseline_3monthsに最尤推定によるガンマ分布を当てはめようとしています。

私がしたいことは、同じことを計算することですが、ラスタスタックを使って計算することです。

私はcalc()機能でこれをやってみました:

f1<-function(x) { library(MASS) return(fitdistr(x,"gamma")) } gamma.parameters<- calc(x = may_baseline_3months,fun = f1)

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function

が、それはうまくいきませんでした。 注:ラスタースタックには4つのレイヤーしかありません。

EDIT

あなたはここにspi

fitdistrが私の主な目標の手順の一部である例のデータをダウンロードすることができます。私は標準降水量指数を計算しようとしています。私はすでに30年間の毎月の降水量の時系列でそれをしました。最後の行は、私は、ラスタ・スタックのための計算に問題を抱えているものであると

data<-read.csv("guatemala_spi.csv",header = T,sep=";") 

dates<-data[,1] 
rain_1month<-data[,2] 
rain_3months<-0 
#Setting the first 2 elements to NA because I'm going to calcule the accumulating the rainfall for 3 month 
for (i in c(1:2)) { 
    rain_3months[i]<-NA 
} 
#Accumulating the rainfall for the rest of the data 
number_of_months<-length(rain_1month) 
for (j in c(3:number_of_months)) 
{ 
rain_3months[j]<-0.0 
for (i in c(0:2)) 
{ 
rain_3months[j] = rain_3months[j] + rain_1month[j-i] 
} 
} 
#Extracting a time-series for the month of interest (May) 
may_rain_3months<-rain_3months[substr(dates,5,6)==”05”] 
dates_may<-dates[substr(dates,5,6)==”05”] 
number_of_years<-length(dates_may) 

#Fitting the gama distribution by maximum likelihood estimation 
start_year<-1971 
end_year<-2010 
start_index<-which(substr(dates_may,1,4)==start_year) 
end_index<-which(substr(dates_may,1,4)==end_year) 
may_baseline_3months<-may_rain_3months[start_index:end_index] 

library(MASS) 
gamma.parameters<-fitdistr(may_baseline_3months[may_baseline_3months>0],"gamma") 

:ここ

私は株式だラインまでの時系列のためのコードです。

ここで私は、ラスタ形式で、これまで持っているものです。

例多層ラスタhere(月額降水量2001 2004年まで、合計で48層)

#Initiating a dates vector 
dates<-c("200101","200102","200103","200104","200105","200106","200107","200108","200109","200110","200111","200112", 
    "200201","200202","200203","200204","200205","200206","200207","200208","200209","200210","200211","200212", 
    "200301","200302","200303","200304","200305","200306","200307","200308","200309","200310","200311","200312", 
    "200401","200402","200403","200404","200405","200406","200407","200408","200409","200410","200411","200412") 
#Initiating a NA raster 
rain_3months_1layer<-raster(nrow=1600, ncol=1673,extent(-118.4539, -34.80395, -50, 30),res=c(0.05,0.05)) 
values(rain_3months_1layer)<-NA 

#Creating a raster stack NA of 48 layers 
rain_3months<-stack(mget(rep("rain_3months_1layer" , 48))) 

#Reading the data 
rain_1month <- stack("chirps_rain_1month.tif") 

#Accumulating the rainfall 
number_of_months<-nlayers(rain_1month) 
for (j in c(3:number_of_months)) 
{ 
rain_3months[[j]]<-0.0 
for (i in c(0:2)) 
{ 
rain_3months[[j]] = rain_3months[[j]] + rain_1month[[j-i]] 
} 
} 

#Extracting the raster for the month of interest (May) 
may_rain_3months<-stack(rain_3months[[which(substr(dates,5,6)=="05", arr.ind = T)]]) 
dates_may<-dates[substr(dates,5,6)=="05"] 
number_of_years<-length(dates_may) 

#Fitting the gama distribution by maximum likelihood estimation 
start_year<-2001 
end_year<-2004 
start_index<-which(substr(dates_may,1,4)==start_year) 
end_index<-which(substr(dates_may,1,4)==end_year) 
may_baseline_3months<-stack(may_rain_3months[[start_index:end_index]]) 

library(MASS) 

f1<-function(x) 
{ 
library(MASS) 
return(fitdistr(x,"gamma")) 
} 
gamma.parameters<- calc(x = may_baseline_3months,fun = f1) 

私はcalc()にすることはできませんラスタスタックにfitdistr()を計算します。

+0

rstudioタグは無関係です(これは普通のRです)、ラスタタグを追加しました。私たちはすべて実行できる再現可能な例を作れますか? – Spacedman

+0

リプレイをありがとう!私は例を追加しました。 –

答えて

0

calcが使用できる機能を作成する必要があります。あなたの関数f1はクラスfitdistrのオブジェクトを返します。 calc機能はそれで何をすべきかを知りません。

library(MASS) 
set.seed(0) 
x <- runif(10) 

f1 <- function(x) { 
    return(fitdistr(x,"gamma")) 
} 

a <- f1(x) 
class(a) 
# [1] "fitdistr" 
a 
# shape  rate 
# 4.401575 6.931571 
# (1.898550) (3.167113) 

あなたは数字を返す関数を必要としています。同様f2calc

f2 <- function(x) { 
    fitdistr(x,"gamma")$estimate 
} 

b <- f2(x) 

class(b) 
#[1] "numeric" 
b 
# shape  rate 
#4.401575 6.931571 

テストf2

library(raster) 
s <- stack(lapply(1:12, function(i) setValues(r, runif(ncell(r))))) 
r <- calc(s, f2) 

私はこれがあなたの質問に答えることを前提としています。あなたの質問があまりにも複雑すぎるので、私は確信できません。このような問題で最初にやるべきことは、私が上で行ったような簡単な例を作ることです。統計に

次の質問

エラー:: OPTIM(図示X = Cの(7、7、7、7)、PAR =リスト(形状= Infを、レート = INF):非OPTIMによって供給-finite値。

別の問題である、あなたはそれに対処することができない値でfitdistrを提供している。あなたがそれらをスキップするtry句を追加することができます。あなたはこれがで起きた細胞を識別でき、あなたが何か他のものがあるなら、値が何を参照するべきかo。あなたはf3で返される値の数は常に同じであることを確認する必要があり

f3 <- function(x) { 
    x <- try (fitdistr(x,"gamma")$estimate, silent=TRUE) 
    if (class(x) == 'try-error') { c(-9999, -9999) } else { x } 
} 

x[1] <- NA 
f2(x) 
#Error in fitdistr(x, "gamma") : 'x' contains missing or infinite values 
f3(x) 
#[1] -9999 -9999 

注意。この場合、2つの値。ここでは、-9999を使用して、細胞を特定することができます。 NA

+0

ありがとう!それはまさに私が欲しいものですが、私はこの問題を抱えています。統計情報:: optim(x = c(7,7,7,7)、par = list(shape = Inf、rate = Inf))のエラー: non-finite最適値によって供給される値。 –

+0

私は上記のあなたのフォローアップ質問への回答を追加しました – RobertH

関連する問題