首页 文章

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

提问于
浏览 次
1

我一直在寻找一种方法来从一系列海面温度（sst）矩阵中按月选择R范围内的值的频率 . 这些数据排列为（sst，80 * 40 * 172），即北大西洋的（sst，经度，纬度，月） . 我对这些频率感兴趣，因为我用它来计算每个月sst范围的表面积，比如在4°C到12°C之间的等温线 . 我经常使用R对时间序列和空间数据进行统计分析，但我不习惯编程，所以我的方法可能不是最有效的 . 我成功地为一个月内的所有纬度提取了> 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）提取，并应用两个自定义函数表面区域的lat条带SA和所有lat条带的表面区域为月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或表面区域是一个纬度（纬度）带SA的表面和月SAM的所有带的面积，用于温度上限（tmu） .

来自函数SA的矩阵的表面积为sst> tmu（比如12°C），纬度（纬度）（比如30°N），而函数SAM的矩阵（sa）具有所有纬度的总和（来自30°N到70°N） . 这可以完成一个月的工作，我可以重复这些功能12次，以16年的方式获得这一年或172次：

我定义了一个起始纬度（纬度），然后是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 . 提前致谢 .

r

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)

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

    回复于

相关问题