R语言气象相关分析方法

R语言气象相关分析方法（包括二进制文件R语言读取）

######读取热带印度洋海温场资料#################
read.filename <- file("D:\\99\\NCEP_IOSST_30y_Wt.dat", "rb")
bindata <- readBin(read.filename, numeric(),16*7*30,size = 4)
IOSST<-array(data = bindata,dim = c(16,7,30))
close(read.filename)
#######读取热带太平洋海温场资料##########
read.filename <- file("D:\\99\\NCEP_TPSST_30y_Wt.dat", "rb")
bindata <- readBin(read.filename, numeric(),32*7*30,size = 4)
TPSST<-array(data = bindata,dim = c(32,7,30))
close(read.filename)
######读取全球200hPa位势高度场资料########
read.filename <- file("D:\\99\\NCEP_Z200_30y_Wt.dat", "rb")
bindata <- readBin(read.filename, numeric(),48*24*30,size = 4)
Z200<-array(data = bindata,dim = c(48,24,30))
close(read.filename)
#######计算IOB海温序列（热带印度洋区域平均海温）#####
IOB<-matrix(0,1,30)
IOSST[IOSST>1e+32]<-NA
for (i in 1:30) {0
  j<-sum(is.na(IOSST[,,i]))
  IOSST[,,i][is.na(IOSST[,,i])]<-0
  IOB[1,i]<-sum(IOSST[,,i])/(16*7-j)
}
#######计算ElNino海温指数######
Nino<-array(TPSST[13:22,3:5,1:30],dim=c(10,3,30))
EY<-matrix(0,1,30)
for (i in 1:30) {
  EY[i]<-sum(Nino[,,i])/30
}
#######一元回归求位势高度场与热带太平洋ElNino海温的关系#####
coeff<-matrix(0,48,24)
Q<-matrix(0,48,24)
c<-matrix(0,48,24)
F<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    relation<-data.frame(X<-as.numeric(EY),Y<-as.numeric(Z200[i,j,]))
    c[i,j]<-1/(sum((relation[,1]-mean(relation[,1]))^2))
    lm.sol<-lm(Y~1+X,data = relation)
    Q[i,j]<-sum((relation[,2]-lm.sol$fitted.values)^2)
    coeff[i,j]<-lm.sol$coefficients[2]
    F[i,j]<-(coeff[i,j]^2/c[i,j])/(Q[i,j]/28)
  }
}
write.filename<-file("D:\\99\\sx3\\coeff_nino.grd","wb")
writeBin(as.vector(coeff),write.filename,size = 4)
close(write.filename)
write.filename<-file("D:\\99\\sx3\\f_nino.grd","wb")
writeBin(as.vector(F),write.filename,size = 4)
close(write.filename)
#########一元回归分析位势高度场与热带印度洋海盆一致模IOB海温异常之间的关系####
coeff<-matrix(0,48,24)
Q<-matrix(0,48,24)
c<-matrix(0,48,24)
F<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    relation<-data.frame(X<-as.numeric(IOB),Y<-as.numeric(Z200[i,j,]))
    c[i,j]<-1/(sum((relation[,1]-mean(relation[,1]))^2))
    lm.sol<-lm(Y~1+X,data = relation)
    Q[i,j]<-sum((relation[,2]-lm.sol$fitted.values)^2)
    coeff[i,j]<-lm.sol$coefficients[2]
    F[i,j]<-(coeff[i,j]^2/c[i,j])/(Q[i,j]/28)
  }
}
write.filename<-file("D:\\99\\sx3\\coeff_IOB.grd","wb")
writeBin(as.vector(coeff),write.filename,size = 4)
close(write.filename)
write.filename<-file("D:\\99\\sx3\\f_IOB.grd","wb")
writeBin(as.vector(F),write.filename,size = 4)
close(write.filename)
#######多元回归分析位势高度场与热带太平洋El Nino海温的关系####
coeff_1<-matrix(0,48,24)
coeff_2<-matrix(0,48,24)
F<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    relation<-data.frame(X1<-as.numeric(EY),X2<-as.numeric(IOB),Y<-as.numeric(Z200[i,j,]))
    lm.sol<-lm(Y~X1+X2,data = relation)
    coeff_1[i,j]<-lm.sol$coefficients[2]
    coeff_2[i,j]<-lm.sol$coefficients[3]
  }
}
write.filename<-file("D:\\99\\sx3\\coeff_1.grd","wb")
writeBin(as.vector(coeff_1),write.filename,size = 4)
close(write.filename)
write.filename<-file("D:\\99\\sx3\\coeff_2.grd","wb")
writeBin(as.vector(coeff_2),write.filename,size = 4)
close(write.filename)
for (j in 1:24) {
  for (i in 1:48) {
    u<-lm.sol$coefficients[2]*cov(as.numeric(EY),as.numeric(Z200[i,j,]))+lm.sol$coefficients[3]*cov(as.numeric(IOB),as.numeric(Z200[i,j,]))
    R_squared<-u/(30*sd(Z200[i,j,]))
    F[i,j]<-(R_squared/2)/((1-R_squared)/(30-2-1))
  }
}
write.filename<-file("D:\\99\\sx3\\dyhg_f.grd","wb")
writeBin(as.vector(F),write.filename,size = 4)
close(write.filename)
######偏相关分析分析位势高度场与热带太平洋El Nino海温的关系####
library(ggm)  
pcor_1<-matrix(0,48,24)
pcor_2<-matrix(0,48,24)
pcor_3<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    relation<-data.frame(X1<-as.numeric(EY),X2<-as.numeric(IOB),Y<-as.numeric(Z200[i,j,]))
    pcor_1[i,j]<-pcor(c(3,2,1),cov(relation))
  }
}
for (j in 1:24) {
  for (i in 1:48) {
    relation<-data.frame(X1<-as.numeric(EY),X2<-as.numeric(IOB),Y<-as.numeric(Z200[i,j,]))
    pcor_2[i,j]<-pcor(c(3,1,2),cov(relation))
  }
}
for (j in 1:24) {
  for (i in 1:48) {
    relation<-data.frame(X1<-as.numeric(EY),X2<-as.numeric(IOB),Y<-as.numeric(Z200[i,j,]))
    pcor_3[i,j]<-pcor(c(1,2,3),cov(relation))
  }
}
write.filename<-file("D:\\99\\sx3\\pcor_1.grd","wb")
writeBin(as.vector(pcor_1),write.filename,size = 4)
close(write.filename)
write.filename<-file("D:\\99\\sx3\\pcor_2.grd","wb")
writeBin(as.vector(pcor_2),write.filename,size = 4)
close(write.filename)
write.filename<-file("D:\\99\\sx3\\pcor_3.grd","wb")
writeBin(as.vector(pcor_3),write.filename,size = 4)
close(write.filename)
#######偏相关分析分析显著性检验,t检验，自由度为27######
pcor_t1<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    a<-pcor.test(pcor_1[i,j],1,30)
    pcor_t1[i,j]<-a$tval
  }
}
write.filename<-file("D:\\99\\sx3\\pcor_t1.grd","wb")
writeBin(as.vector(pcor_t1),write.filename,size = 4)
close(write.filename)
###
pcor_t2<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    a<-pcor.test(pcor_2[i,j],1,30)
    pcor_t2[i,j]<-a$tval
  }
}
write.filename<-file("D:\\99\\sx3\\pcor_t2.grd","wb")
writeBin(as.vector(pcor_t2),write.filename,size = 4)
close(write.filename)
###
pcor_t3<-matrix(0,48,24)
for (j in 1:24) {
  for (i in 1:48) {
    a<-pcor.test(pcor_3[i,j],1,30)
    pcor_t3[i,j]<-a$tval
  }
}
write.filename<-file("D:\\99\\sx3\\pcor_t3.grd","wb")
writeBin(as.vector(pcor_t3),write.filename,size = 4)
close(write.filename)