 登録日: 2004-7-29 居住地: 地球 投稿: 303 |
定義に沿ったAOの計算 北極振動を TW98 に沿って計算し直す。 1. CD-ROM データからの切り出し #! /bin/bash
ROOTDIR=/DATA/DL/NCEP_NCAR/
FILENAM=data/monthly/other/pres.msl
#
for dir in ${ROOTDIR}??
do
wgrib $dir/$FILENAM | grep 'd=..[01][1-4]' | wgrib $dir/$FILENAM -i -bin -nh -o ${dir##*/}.bin
done
各年 252288 = 144 x 73 x 6 x 4 バイトのファイルができる。 2. R で読み込んで保存する。 2.1 取り込み # 1. データの読み込み
#
data <- NULL
for ( year in c(50:99, "00", "01", "02" )){
zz <- file(paste(year,".bin",sep=""), "rb")
data <- c( data, readBin(zz, "double", n=3564928, size=4)) # 多めに指定
close(zz)
}
rm(zz)
dim(data) <- c(144,73,6,53)
data <- data[,seq(73,1,-1),,]
data <- data[,seq(45,73),,]
2.2 とりあえずテスト表示 source("~/lib/R/MAP/map.R")
source("~/lib/R/filledcontour.R")
filled.contour(seq(0,357.5,2.5),seq(20,90,2.5),data[,,1,1])
map("/home/mori/lib/R/MAP/coast_world.asc")
2.3 保存 save(data, file="NHSLPWIN.rdata")
3. 解析 3.1 アノマリ
for( month in 1:6 ){
meandata <- data[,,1,1]*0
for( year in 1:53 ){
meandata <- meandata + data[,,month,year]
}
meandata <- meandata / 53
for( year in 1:53 ){
data[,,month,year] <- data[,,month,year] - meandata
}
}
3.2 重みづけ data <- data[,seq(1,28),,]
weight <- diag(sqrt(cos( ( 20 + seq(0,27)*2.5 ) /180*pi ) ))
for( month in 1:6 ){
for( year in 1:53 ){
data[,,month,year] <- data[,,month,year] %*% weight
}
}
3.3 主成分分析 dim(data) <- c(144*28,6*53)
EOFresult <- prcomp(t(data))
save(EOFresult, file="EOFresult.rdata")
3.4 主成分分析の結果の吟味 #
# http://jisao.washington.edu/analyses0302/ との比較
#
refdata <- t(matrix(scan("slpanompc.ascii.txt"),5))
AOI <- EOFresult$x[,1] ; dim(AOI) <- c(6,53)
PC2 <- EOFresult$x[,2] ; dim(PC2) <- c(6,53)
PC3 <- EOFresult$x[,3] ; dim(PC3) <- c(6,53)
plot(AOI[1,], refdata[refdata[,2]==1,][3:55,3])
plot(AOI[2,], refdata[refdata[,2]==2,][3:55,3])
plot(PC2[1,], refdata[refdata[,2]==1,][3:55,4])
plot(PC3[1,], refdata[refdata[,2]==1,][3:55,5])
※ 今回の場合、第一成分の符号を変えるべきだということがわかる。
# PC1 について符号を変更
EOFresult$x[,1] <- - EOFresult$x[,1]
EOFresult$rot[,1] <- - EOFresult$rot[,1]
# Pa -> hPa
EOFresult$x <- EOFresult$x/100
3.5 空間分布の表示 source("00_Functions.R")
eigenplot(EOFresult$rot[,1])
eigenplot(EOFresult$rot[,2])
ここで、00_Functions.R の中で、次のような関数を定義している。 ###################################################################################
#
# Polar の等値線図を描く関数
#
eigenplot <- function( eigenvector, zlim=c(-0.1,0.1), title="" )
{
library(akima)
phi <- matrix( seq( 0,357.5, 2.5),144,28)
the <- t(matrix(seq( 20, 87.5, 2.5),28,144))
weight <- diag(sqrt( cos( 2 * pi / 360 * ( 20 + seq(0,27)*2.5 )) ))
polardata <- interp((90-the) * ( -sin(2*pi/360 * phi) ),
(90-the) * ( cos(2*pi/360 * phi) ),
matrix(eigenvector,144) %*% diag(1./diag(weight)))
filled.contour(polardata, zlim=zlim, nlevels=15, plot.axes={} )
contour( polardata, zlim=zlim, nlevels=15, add=T )
###############
if( is.na(file.info("/home/mori/lib/R/MAP/coast_world.rdata")$size) )
saveworldmapdata()
load("/home/mori/lib/R/MAP/coast_world.rdata")
###############
for( i in seq(1,mapdata$max) )
lines( -(90 - mapdata$data[[i]]$lat) * sin(2*pi*mapdata$data[[i]]$lon/360),
( 90 - mapdata$data[[i]]$lat) * cos(2*pi*mapdata$data[[i]]$lon/360))
title(title)
}
###################################################################################
#
# 汎用
# 世界地図の海岸線データを保存するための関数
#
saveworldmapdata <- function(file="/home/mori/lib/R/MAP/coast_world.asc"){
buff <- scan(file)
max <- length(buff)
sfx <- 1
cnt <- 0
mapdata <- list(list())
while( sfx <= max ){
cnt <- cnt + 1
length <- buff[sfx]
mapdatatmp <- list(num = buff[sfx]/2,
unknown = buff[sfx+1],
north = buff[sfx+2],
south = buff[sfx+3],
east = buff[sfx+4],
west = buff[sfx+5],
lat = buff[sfx+4+2*c(1:(length/2))],
lon = buff[sfx+5+2*c(1:(length/2))])
mapdata[[cnt]] <- mapdatatmp
sfx <- sfx + 6 + length
}
mapdata <- list(max=cnt,data=mapdata)
save(mapdata, file=("/home/mori/lib/R/MAP/coast_world.rdata"))
}
4. 独立成分分析 4.1 やるだけ。 library(fastICA)
ICAresult <- fastICA(EOFresult$x[,1:3],3,v=T,tol=1e-14,maxit=2000)
4.2 空間分布の表示 eigenplot( (EOFresult$rot[,1:3] %*% t(ICAresult$A))[,1],zlim=c(-1000,1000) )
規格化などをしたので、独立成分1σに対する気圧変動量(1hPa)の図になる。 |
メイン