1

我一直在寻找一种方法来按月从一系列海面温度 (sst) 矩阵中选择 R 范围内的值的频率。这些数据排列为 (sst, 80*40*172)，即北大西洋的 (sst, longitude,latitude,month)。我对这些频率很感兴趣，因为我使用它们来计算每个月 sst 范围内的表面积，例如 4°C 和 12°C 等温线之间。我经常使用 R 对时间序列和空间数据进行统计分析，但我不习惯编程，所以我的方法可能不是最有效的。我已经成功地提取了一个月内所有纬度的 sst 值 > 12 °C 的频率：

我来自 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)，并应用两个自定义函数表面区域的纬度带 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 或表面积 一个纬度 (lat) 条带 SA 的表面和月份 SAM 的所有条带的面积，用于温度上限 (tmu)。

来自函数 SA 的矩阵 (s) 对于纬度 (lat)（比如在 30°N）具有 sst> tmu（比如 12 °C）的表面积，来自函数 SAM 的矩阵 (sa) 具有所有纬度的总和（从30°N 到 70°N）。这完成了一个月的工作，我可以重复这些功能 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 回答