sst数据矩阵序列的循环或函数

Question

我一直在寻找一种方法，以按月从一系列海面温度 (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) 的表面积为 sst> tmu(比如 12 °C) 用于纬度 (lat)（比如 30°N），来自函数 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。提前致谢。

Answer 1

如果我理解您的问题，您的代码会简化为这样

#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)

如果这不正确，请让您的问题更清楚，并提供reproducible 输入和输出数据集。