1

月ごとの海面水温（sst）の行列から、Rの範囲内の値の頻度を選択する方法を探していました。これらのデータは、北大西洋の（sst、80 * 40 * 172）、つまり（sst、経度、緯度、月）として配置されます。私はこれらの周波数に興味があります。なぜなら、それらを使用して、4°Cから12°Cの等温線などの毎月のsstの範囲の表面積を計算するからです。私は時系列および空間データの統計分析にRを頻繁に使用しましたが、プログラミングに慣れていないため、私のアプローチはおそらく最も効率的ではありません。私は、1か月ですべての緯度で>12°Cのsst値の周波数を抽出することに成功しました。

IRI/LDEO気候研究ウェブサイトからの私のデータ

dat <- read.table("c:/Temp/[Y+X]datatable.tsv", header=FALSE, sep="\t")

最初の月の5x5サンプルマトリックス

    Nov 1981    30.5000 31.5000 32.5000 33.5000
    9.50000     21.7906 21.9431 22.1324 22.1662
    8.50000     21.7267 21.8573 21.9981 21.8757
    7.50000     21.6644 21.7781 21.8960 21.7393
    6.50000     21.5989 21.7025 21.8044 21.6304

基本的に経度をスキップし、行ラベルも後で日付を追加する予定です。最初に、最初の月を行列（t）として抽出し、2つのカスタム関数の緯度ストリップSAの表面積と月MSAのすべての緯度ストリップの表面積を適用します。

ro=2:81
co=2:41
t <- as.matrix(dat[ro,co])

それで

SA <- function (lat,tmu){
    l <- c(t[,lat]>=tmu,0)
    la <- as.data.frame(l)
    x  <- la[,1] 
    n  <- length(x)
    sau <- array(0,n)
    x. <- lat
    for (i in 1:n) sau[i] <- (x[i]*111.320)*(cos(x.*(3.1415/180))*111.320)
    s <- as.matrix(sum(sau))
}
MSA <- function(tmu){
    m <-1:40
    su <- array(0,0)
    for (i in 1:40) su[i] <-SA(m[i],tmu)
    ms <- as.data.frame(su)
    sa <- as.data.frame(colSums (ms))
    return(sa)
}

関数SAおよびMSAまたは表面積1緯度（緯度）ストリップSAの表面、および温度上限（tmu）の月SAMのすべてのストリップの面積。

関数SAの行列（s）は、緯度（lat）（たとえば30°N）に対してsst> tmu（たとえば12°C）の表面積を持ち、関数SAMの行列（sa）は、すべての緯度の合計（ 30°Nから70°N）これは1か月間機能し、関数を12回繰り返して年を取得するか、次のように16年間で172回繰り返すことができます。

開始緯度（lat）を定義し、次に81経度セル（st）のステップを定義して、次に抽出し、目的の温度を定義します。

lat= 30.5 
st=81 
t=12

次に、各月の総表面積を計算します。

SA1 <- {
    i=0*st 
    t <- as.matrix(dat[ro+i,co])
    SA(lat,t)
    MSA(12)
}

SA2 <- {
    i=1*st 
    t <- as.matrix(dat[ro+i,co])
    SA(30.5,t)
    MSA(12)
}

私の質問は、172回のすべての月を反復するループまたは関数を作成できるかどうかです。したがって、SA1、SA2...SA172の繰り返しをスキップします。前もって感謝します。

4

1 に答える 1

2

私があなたの質問を理解すれば、あなたのコードはこれに還元されます

#create some dummy data
lat <- 30:33
lon <- 6:9
sst <- array(rnorm(length(lat) * length(lon) * 12, mean = 21, sd = 4), 
    dim = c(length(lat), length(lon), 12), dimnames = list(lat, lon, 1:12))
#the actual code
SA <- function (test, tmu){
  colSums(test >= tmu) * 111.320 ^ 2 * 
     cos(as.numeric(rownames(test)) * (pi/180))
}
apply(sst, 3, SA, tmu = 21)

これが正しくない場合は、質問をより明確にし、再現可能な入力および出力データセットを提供してください。

于 2012-07-25T13:33:19.470 に答える