argo-sine-fitting

Here are a few plots derived from APDRC ARGO data.  They were obtained by fitting a sine curve to 12 months of potential temperature data to obtain the phase, amplitude, and mean of the annual signal.  Climatology data is plotted in addition to the six years of available ARGO data.

The sampled data was taken from the latitudinal slice at 45ºS.  I chose this latitude as I surmised that it had the advantages of a sinusoidal solar forcing, good amplitude, a high ocean to land ratio, and it being midway between the skewing effects of the tropics and poles.  From a Newtonian view, phase shift and amplitude are measures of how quickly temperatures equalize.  Possibly this might be a quick and dirty way of assessing the heat uptake into the oceans.  Note that I have no experience in thermal dynamics or oceanography, so feel free to disregard anything you see here.

First up... phase relative to the September equinox.  Based on climatology temperatures lag the solar forcing by about 70 days at the surface.  The phase then extends out to about 150 days by 250 M depth, then remains fairly constant to 1000 M although there is a slight reduction.  Not surprisingly, the climatology line in much more regular than those of the individual years.  There is a lot of variation in the individual years and an interesting oscillation between a short phase one year and a long phase the next.  Will this pattern persist?  I can't think of a good reason for it, so I doubt it.  From a coin flipping analogy, however, there is a 1 in 32 chance of getting this pattern in six flips.

Next up... potential temperature amplitude.  The graph is showing about a 1/4 degree variation at the surface.  More interestingly, after 75M it becomes apparent that the long phase years have a lesser amplitude relative to climatology and the short phase years a slightly higher amplitude.  This is consistent with the Newtonian view.

Next there's mean temperature.  Not much to see here.  Possible you can see the mean trending upwards, but this plot isn't well suited for trend analysis.

Here's the Adjusted R^2 verification statistic for the linear regression.  For the individual years this is ugly, with some years getting into negative territory at a fairly shallow depth.  Based upon this, any observation I made earlier regarding individual years could be easily dismissed.  Perhaps the results were due to the loss of signal to noise?  The climatology line maintains a reasonable score, so it might be suitable to compare against the climate models.

Here's a plot of the fitted curve to the actual signal for Climatology data.  The APDRC data is derived using a process called "variational interpolation".  Could the amplitude signal be partially transferred between the different layers?  The pattern between 400M and 1000M looks amazingly similar.  Note the wedge in the top of the wave between 250M and 1000M.  Will this shape persist, diminish, or grow?  If it doesn't persist, what effect will it have on the amplitude/signal?  If it does persist, what causes it?  Is there a annual/bi-annual oscillation of upwelling cold water from the deep?

For comparison purposes, here's a plot of the fitted curve to the actual signal for year 2008 data.  Does this model underestimate the amplitude once the signal starts to become lost?


Here's the results of fitting a sine curve at all latitudinal slices between 55º North and South.  The longitudinal dimension was collapsed by taking the mean of all points for each latitude.  For the Lag Plot, by visual inspection, I switched from a September to a March equinox offset at 7ºN.  There is a half year insolation period at the equator, so an annual fit is not expected in this region.  




Here's the R code for the programs that produced the above graphs:

#
# ARGO Sine Fitting
#
# This program fits a sine curve to 12 months of potential temperature data at different depths
#
# The programs takes a latitudinal slice of data from 45ºS and averages the temperature by month and depth.
# A sine curve is then fitted at each depth to get the phase, amplitude, and mean for the annual signal.
# Graphs of these values and the R^2 verification statistic are then plotted.
# Two netCDF files of ARGO data are used from the Asia-Pacific Data-Research Center (apdrc)
#  1x1 gridded Monthly climatology netCDF Data file:
#   http://apdrc.soest.hawaii.edu/projects/Argo/data/gridded/On_standard_levels/argo_CLIM_grd.nc
#  2005-2012 1x1 gridded Monthly means netCDF Data file:
#   http://apdrc.soest.hawaii.edu/projects/Argo/data/gridded/On_standard_levels/argo_2005-2012_grd.nc
#
# Note: I'm new to R, so my apologies in advance.
#

require(ncdf)
require(abind)

nc1<-open.ncdf("argo_CLIM_grd.nc")
nc2<-open.ncdf("argo_2005-2012_grd.nc")

#
#### PARAMETERS ####

sy<-01 # 1st year = 1
cy<-06 # number of years

ub<-0      # Upper Level Depth
lb<-1050   # Lower Level Depth - ARGO data goes to 2000M

wb<-0       # western boundary
eb<-360     # eastern boundary

samplat<- -45

nb<- samplat + 0.5
sb<- samplat - 0.5

#
# Determine the phase offset days for the desired equinox.
# My calculations may be off by a day or two.
#

hem<-"S"

if (hem=="N") {pod<- -80.25} else {pod<- 102.25}

#
# Function: getncsample
# This function returns a four dimensional array of netCDF potential temperature data.
# In addition to the nc file and source description, the function also
# receives the bounds for the sample:
#  wb: west bound
#  eb: east bound
#  sb: south bound
#  nb: north bound
#  ub: upper bound (ocean depth - usually zero)
#  lb: lower bound
#  sy: start year
#  cy: count years (number of years to sample)
#

getncsample<-function(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy)
{
 if (wb <= eb){
    ptsamp = selectncsamp(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy)
 } else 
 {
    ptsamp1 = selectncsamp(nc, ncsource, wb, 360, sb, nb, ub, lb, sy, cy)
 ptsamp2 = selectncsamp(nc, ncsource,  0,  eb, sb, nb, ub, lb, sy, cy)
    ptsamp  = abind(ptsamp1, ptsamp2, along=1)
 }
 return(ptsamp)
}

#
# Function: selectncsample
# Called from getncsample, this does the actual data retrieval from
# the nc file.  Same parameters as getncsample
#

selectncsamp <- function(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy){

 if (ncsource=="argo"){ 
    longvar<-"LONGITUDE"
    latvar<-"LATITUDE"
    depvar<-"LEVEL"
    timevar<-"TIME"
    tempvar<-"PTEMP"
    k<-0} else
    if (ncsource=="standard1"){ 
       longvar<-"lon"
       latvar<-"lat"
       depvar<-"depth"
       timevar<-"time"
       tempvar<-"thetao"
       k<-273.15}
 
 long<-get.var.ncdf(nc,longvar)
 lat<-get.var.ncdf(nc,latvar)
 depth<-get.var.ncdf(nc,depvar)
 
 longx<-which(long>=wb & long<=eb)
 latx<-which(lat>=sb & lat<=nb)
 depthx<-which(depth>=ub & depth<=lb)
 
 startlong<-min(longx)
 countlong<-length(longx)
 
 startlat<-min(latx)
 countlat<-length(latx)
  
 startdepth<-min(depthx)
 countdepth<-length(depthx)
 
 sm=(sy-1)*12+1
 cm=cy*12
 
 start<-c(startlong,startlat,startdepth,sm)
 count<-c(countlong,countlat,countdepth,cm)
 
 ptsamp<-get.var.ncdf(nc,tempvar,start,count) - k
 dim(ptsamp)<-count

 attributes(ptsamp)$dimnames <- list(long[longx],lat[latx],depth[depthx],1:cm-0.5)

 return(ptsamp)
 
}

#
# Function: sinefit
# This function fits a sine wave to 12 months of data and returns
# the mean, phase, amplitude, and R^2 verification statistic.
#
# The input variable y is a two dimensional array with rows representing depths
# and columns representing months.  The variable phaseoffsetdays is used to mark
# the desired equinox.
#

sinefit<-function(y,phaseoffsetdays)
{
 y<-t(y)
 x<-(2*pi*(((0:11+.5)/12)+(phaseoffsetdays/365)))
 
 x1<-cos(x)
 x2<-sin(x)
 
 m<-lm(y ~ x1 + x2)
 s<-summary(m)
 r<-coef(m)
 a<-r[1,]
 b<- -atan2(r[2,],r[3,])
 b<-ifelse(b<0,2*pi+b,b)
 b<-(b*365)/(2*pi)
 c<-sqrt(r[2,]^2+r[3,]^2)
 ############ adjust phase
 
 bdiff = diff(b)
 
 x = bdiff < -(365/2)
 bdiff[x] = bdiff[x] + 365
 x = bdiff >= (365/2)
 bdiff[x] = bdiff[x] - 365
 
 bdiff = c(b[1],bdiff)
 
 b = cumsum(bdiff)
 
 ############
 wmean<-a
 wphase<-b
 wamp<-c
 wr2<-NULL
 for (i in 1:length(s)) {wr2<-c(wr2,s[[i]]$adj.r.squared)} 
 ret<-cbind(wmean,wphase,wamp,wr2)
 return(ret)
}

#
#### MAIN PROGRAM ####

#
# Get the sample of climatology data (sy=1, cy=1).
#

ncsource1<-"argo"
ptnc1<-getncsample(nc1, ncsource1, wb, eb, sb, nb, ub, lb, 1, 1)  

#
# Get rid of longitudes that have NA values.
#

ptlong1<-apply(ptnc1,1,function(x)all(!is.na(x)))
ptlong1<-which(ptlong1==TRUE)
ptsamp1<-ptnc1[ptlong1,,,]

#
# Calculate averages for depth and month, then fit sine.
#

ptmean1<-apply(ptsamp1,c(3,4),mean)
ptfit1<-sinefit(ptmean1,pod)

#
# Save depths for plotting.
#

ptdepth1<-as.numeric(dimnames(ptsamp1)[[3]])

#
# Get the sample of yearly data.
#

ncsource2<-"argo"
ptnc2<-getncsample(nc2, ncsource2, wb, eb, sb, nb, ub, lb, sy, cy)  

#
# Reorganize month dimension into month and year dimensions.
#

d<-dim(ptnc2)
dn<-dimnames(ptnc2)
dim(ptnc2)<-c(d[1],d[2],d[3],12,cy)
attributes(ptnc2)$dimnames <- list(dn[1][[1]],dn[2][[1]],dn[3][[1]],0:11+0.5,1:cy)

#
# Get rid of longitudes that have NA values.
#

ptlong2<-apply(ptnc2,1,function(x)all(!is.na(x)))
ptlong2<-which(ptlong2==TRUE)
ptsamp2<-ptnc2[ptlong2,,,,]

#
# Save depths for plotting.
#

ptdepth2<-as.numeric(dimnames(ptsamp2)[[3]])

#
# Start plotting four graphs
#

dev.new()
dphase<-dev.cur()
dev.new()
damp<-dev.cur()
dev.new()
dmean<-dev.cur()
dev.new()
dr2<-dev.cur()

#
# Plot climatology data
#

###### PHASE ########

dev.set(dphase)

xmin<- 0
xmax<- 400
ymin<-floor(min(ptdepth1))
ymax<-ceiling(max(ptdepth1))

plot(NA,NA,type="n",main=c("ARGO Phase Relative to Sept. Equinox at 45ºS","Climatology vs Individual Yrs"),
ylab="Depth (M)", xlab="Phase Difference (Days)",
  xlim=c(xmin, xmax), ylim=c(ymax, ymin))

abline(h=seq(0,1000,50),v=seq(0,375,25),lty="dotted",col="lightgrey")

lines(ptfit1[,"wphase"], ptdepth1, type = "o", pch=19, col=1, lwd=2)

#
# I will be creating a gif of phase plots, so each frame needs to be saved.
#

#filename<-"argo-phase-45s-0.jpg"
#dev.copy(jpeg,filename,width = 450, height = 450, quality=100)
#dev.off()
 
###### AMPLITUDE ########

dev.set(damp)

xmin<-floor(min(ptdepth1))
xmax<-ceiling(max(ptdepth1))
ymin<-floor(min(ptfit1[,"wamp"]))
ymax<-ceiling(max(ptfit1[,"wamp"]))

plot(NA,NA,type="n",main=c("ARGO PTemp Amplitude at 45ºS","Climatology vs Individual Yrs"),
xlab="Depth (M)", ylab="Amplitude (ºC)",
  xlim=c(xmin, xmax), ylim=c(ymin, ymax))

abline(h=seq(0,ymax,0.25),v=seq(0,1000,50),lty="dotted",col="lightgrey")

lines(ptdepth1, ptfit1[,"wamp"], type = "o", pch=19, col=1, lwd=2)

###### MEAN ########

dev.set(dmean)

xmin<-floor(min(ptdepth1))
xmax<-ceiling(max(ptdepth1))
ymin<-floor(min(ptfit1[,"wmean"]))
ymax<-ceiling(max(ptfit1[,"wmean"]))

plot(NA,NA,type="n",main=c("ARGO PTemp Mean at 45ºS","Climatology vs Individual Yrs"),
xlab="Depth (M)", ylab="Mean (ºC)",
  xlim=c(xmin, xmax), ylim=c(ymin, ymax))

abline(h=seq(0,ymax,0.5),v=seq(0,1000,50),lty="dotted",col="lightgrey")

lines(ptdepth1, ptfit1[,"wmean"], type = "o", pch=19, col=1, lwd=2)

###### R2 ########

dev.set(dr2)

ymin<- -0.5
ymax<- 1
xmin<-floor(min(ptdepth1))
xmax<-ceiling(max(ptdepth1))

plot(NA,NA,type="n",main=c("ARGO R2 Statistic For Regression at 45ºS","Climatology vs Individual Yrs"),
xlab="Depth (M)", ylab="Adjusted R Squared",
  xlim=c(xmin, xmax), ylim=c(ymin, ymax))

abline(h=seq(ymin,ymax,0.1),v=seq(0,1000,50),lty="dotted",col="lightgrey")

lines(ptdepth1, ptfit1[,"wr2"], type = "o", pch=19, col=1, lwd=2)

#
# Plot yearly data.  Loop through the years.
#

for (i in 1:cy)
{

#
# Calculate averages for depth and month, then fit sine.
#
  ptmean2<-apply(ptsamp2[,,,,i],c(3,4),mean)
  ptfit2<-sinefit(ptmean2,pod)

###### PHASE ########

  dev.set(dphase)
  
  lines(ptfit2[,"wphase"], ptdepth2, type = "o", pch=i+1, col=i+1)

#
# Save phase frame for gif.
#
#  filename<-paste("argo-phase-45s-",i,".jpg",sep="")
#  dev.copy(jpeg,filename,width = 450, height = 450, quality=100)
#  dev.off()

###### AMPLITUDE ########

  dev.set(damp)
  
  lines(ptdepth2, ptfit2[,"wamp"], type = "p", pch=i+1, col=i+1)

###### MEAN ########

  dev.set(dmean)
  
  lines(ptdepth2, ptfit2[,"wmean"], type = "o", pch=i+1, col=i+1)

###### R2 ########

  dev.set(dr2)
  
  lines(ptdepth2, ptfit2[,"wr2"], type = "o", pch=i+1, col=i+1)

}

#
# Add legend and save plots.
#

dev.set(dphase)

legend("topright",inset=0.01,legend=c("CLIM",(2004+sy):(2004+cy)),bg="white",
         text.col=1,pch=c(19,2:(cy+1)),col=c(1:(cy+1)),lwd=c(2,rep(1,cy)))

filename<-"argo-phase-45s.jpg"
dev.copy(jpeg,filename,width = 450, height = 450, quality=100)
dev.off()

dev.set(damp)

legend("topright",inset=0.01,legend=c("CLIM",(2004+sy):(2004+cy)),bg="white",
         text.col=1,pch=c(19,2:(cy+1)),col=c(1:(cy+1)),lwd=c(2,rep(1,cy)))

filename<-"argo-amp-45s.jpg"
dev.copy(jpeg,filename,width = 450, height = 450, quality=100)
dev.off()

dev.set(dmean)

legend("topright",inset=0.01,legend=c("CLIM",(2004+sy):(2004+cy)),bg="white",
         text.col=1,pch=c(19,2:(cy+1)),col=c(1:(cy+1)),lwd=c(2,rep(1,cy)))

filename<-"argo-mean-45s.jpg"
dev.copy(jpeg,filename,width = 450, height = 450, quality=100)
dev.off()

dev.set(dr2)

legend("bottomleft",inset=0.01,legend=c("CLIM",(2004+sy):(2004+cy)),bg="white",
         text.col=1,pch=c(19,2:(cy+1)),col=c(1:(cy+1)),lwd=c(2,rep(1,cy)))

filename<-"argo-r2-45s.jpg"
dev.copy(jpeg,filename,width = 450, height = 450, quality=100)
dev.off()

close.ncdf(nc1)
close.ncdf(nc2)

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

#
# ARGO Plots Sine Curves
#

require(ncdf)
require(abind)

#
#### PARAMETERS ####

#nc1<-open.ncdf("argo_2005-2012_grd.nc") # uncomment and set sy=04 for 2008
nc1<-open.ncdf("argo_CLIM_grd.nc")

sy<-01 # 1st year = 1
cy<-01 # number of years

ub<-0      # Upper Level Depth
lb<-1050   # Lower Level Depth - ARGO data goes to 2000M

wb<-0       # western boundary
eb<-360     # eastern boundary
hem<-"S"    # hemisphere="N","S"

samplat<- -42.5

nb<- samplat + 0
sb<- samplat - 0

#
# Determine the phase offset days for the desired equinox.
# My calculations may be off by a day or two.
#

if (hem=="N") {pod<- -80.25} else {pod<- 102.25}

#
# Function: getncsample
# This function returns a four dimensional array of netCDF potential temperature data.
# In addition to the nc file and source description, the function also
# receives the bounds for the sample:
#  wb: west bound
#  eb: east bound
#  sb: south bound
#  nb: north bound
#  ub: upper bound (ocean depth - usually zero)
#  lb: lower bound
#  sy: start year
#  cy: count years (number of years to sample)
#

getncsample<-function(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy)
{
 if (wb <= eb){
    ptsamp = selectncsamp(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy)
 } else 
 {
    ptsamp1 = selectncsamp(nc, ncsource, wb, 360, sb, nb, ub, lb, sy, cy)
 ptsamp2 = selectncsamp(nc, ncsource,  0,  eb, sb, nb, ub, lb, sy, cy)
    ptsamp  = abind(ptsamp1, ptsamp2, along=1)
 }
 return(ptsamp)
}

#
# Function: selectncsample
# Called from getncsample, this does the actual data retrieval from
# the nc file.  Same parameters as getncsample
#

selectncsamp <- function(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy){

 if (ncsource=="argo"){ 
    longvar<-"LONGITUDE"
    latvar<-"LATITUDE"
    depvar<-"LEVEL"
    timevar<-"TIME"
    tempvar<-"PTEMP"
    k<-0} else
    if (ncsource=="standard1"){ 
       longvar<-"lon"
       latvar<-"lat"
       depvar<-"depth"
       timevar<-"time"
       tempvar<-"thetao"
       k<-273.15}
 
 long<-get.var.ncdf(nc,longvar)
 lat<-get.var.ncdf(nc,latvar)
 depth<-get.var.ncdf(nc,depvar)
 
 longx<-which(long>=wb & long<=eb)
 latx<-which(lat>=sb & lat<=nb)
 depthx<-which(depth>=ub & depth<=lb)
 
 startlong<-min(longx)
 countlong<-length(longx)
 
 startlat<-min(latx)
 countlat<-length(latx)
  
 startdepth<-min(depthx)
 countdepth<-length(depthx)
 
 sm=(sy-1)*12+1
 cm=cy*12
 
 start<-c(startlong,startlat,startdepth,sm)
 count<-c(countlong,countlat,countdepth,cm)
 
 ptsamp<-get.var.ncdf(nc,tempvar,start,count) - k
 dim(ptsamp)<-count

 attributes(ptsamp)$dimnames <- list(long[longx],lat[latx],depth[depthx],1:cm-0.5)

 return(ptsamp)
 
}

#
# Function: sinefit
# This function fits a sine wave to 12 months of data and returns
# the mean, phase, amplitude, and R^2 verification statistic.
#
# The input variable y is a two dimensional array with rows representing depths
# and columns representing months.  The variable phaseoffsetdays is used to mark
# the desired equinox.
#

sinefit<-function(y,phaseoffsetdays)
{
 y<-t(y)
 x<-(2*pi*(((0:11+.5)/12)+(phaseoffsetdays/365)))
 
 x1<-cos(x)
 x2<-sin(x)
 
 m<-lm(y ~ x1 + x2)
 s<-summary(m)
 r<-coef(m)
 a<-r[1,]
 b<- -atan2(r[2,],r[3,])
 b<-ifelse(b<0,2*pi+b,b)
 b<-(b*365)/(2*pi)
 c<-sqrt(r[2,]^2+r[3,]^2)
 wmean<-a
 wphase<-b
 wamp<-c
 wr2<-NULL
 for (i in 1:length(s)) {wr2<-c(wr2,s[[i]]$adj.r.squared)} 
 ret<-cbind(wmean,wphase,wamp,wr2)
 return(ret)
}

#
#### MAIN PROGRAM ####

#
# Get the sample of data
#

ncsource1<-"argo"
ptnc1<-getncsample(nc1, ncsource1, wb, eb, sb, nb, ub, lb, sy, cy)  

#
# Get rid of longitudes that have NA values.
#

ptlong1<-apply(ptnc1,1,function(x)all(!is.na(x)))
ptlong1<-which(ptlong1==TRUE)
ptsamp1<-ptnc1[ptlong1,,,,drop=FALSE]

#
# Calculate averages for depth and month, then fit sine.
#

ptmean1<-apply(ptsamp1,c(3,4),mean)
ptfit1<-sinefit(ptmean1,pod)

#
# Save depths for plotting.
#

ptdepth1<-as.numeric(dimnames(ptsamp1)[[3]])

#

m<-(2*pi*(((0:11+.5)/12)+(pod/365)))
m<-rep(m,length(ptdepth1))
dim(m)<-c(12,length(ptdepth1))
m=t(m)

f<-ptfit1[,"wmean"]+ptfit1[,"wamp"]*sin((-2*pi*ptfit1[,"wphase"]/365)+m)
dimnames(f)<-dimnames(ptmean1)

x<-0:23+.5
r<-dim(f)[1]
d<-dimnames(f)[[1]]
dev.new(width=10,height=5)
par(mfrow=c(ceiling(r/5),5)) # all plots on one page 
for(i in 1:r){ 
  heading = paste(paste("Depth=",d[i],"M",sep=""), paste("R2=",round(ptfit1[i,"wr2"],2),sep=""),sep="  ")

  yf<-c(f[i,],f[i,]) 
  ya<-c(ptmean1[i,],ptmean1[i,]) 

  ymin<-min(ya,yf)
  ymax<-max(ya,yf)
 
  y<-yf
  plot(x, y, type="n", main=heading, ylim=c(ymin, ymax), ylab="Ptemp (ºC)", xlab="Month")
  lines(x,y,type="l",col=1) 

  y<-ya
  lines(x,y,type="l",col=2) 
}

close.ncdf(nc1)

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

#
# Model Sine Fit Analysis - mean latitude then fit
#

library("RNetCDF")
library("abind")

#
#### PARAMETERS ####

outfile<-"pt-55s-to-55n-mtf.RData"

sy<-01 # 1st year = 1

ub<-0      # Upper Level Depth
lb<-1100   # Lower Level Depth - ARGO data goes to 2000M

wb<-0      # western boundary
eb<-360    # eastern boundary

samplat<- 0

nb<- samplat + 55
sb<- samplat - 55

#ncfiles<-c(
#   "argo_CLIM_grd.nc"),
#   "pcmdi.ipcc4.cccma_cgcm3_1_t63.sresa1b.run1.monthly.thetao_o1_sresa1b_1_cgcm3.1_t63_2001_2025.nc",
#   "pcmdi.ipcc4.gfdl_cm2_1.sresa1b.run1.monthly.thetao_O1.200101-201012.nc",
#   "pcmdi.ipcc4.giss_model_e_r.sresa1b.run1.monthly.thetao_O1.GISS1.SRESA1B.run1.nc",
#   "pcmdi.ipcc4.iap_fgoals1_0_g.sresa1b.run1.monthly.thetao_O1_2000_2009.nc",
#   "pcmdi.ipcc4.inmcm3_0.sresa1b.run1.monthly.thetao_O1_1.nc",
#   "pcmdi.ipcc4.miroc3_2_hires.sresa1b.run1.monthly.thetao_O1_y2001-y2010.nc",
#   "pcmdi.ipcc4.mpi_echam5.sresa1b.run1.monthly.thetao_O1_2001-2009.nc",
#   "pcmdi.ipcc4.ncar_ccsm3_0.sresa1b.run2.monthly.thetao_O1.SRESA1B_2.CCSM.ocnm.2005-01_cat_2009-12.nc",
#   "pcmdi.ipcc4.ncar_pcm1.sresa1b.run2.monthly.thetao_O1.SRESA1B_2.PCM1.ocnm.2000-01_cat_2009-12.nc",
#   "pcmdi.ipcc4.ukmo_hadcm3.sresa1b.run1.monthly.thetao_O1_2000_to_2024.nc"
#)

ncfiles<-"argo_CLIM_grd.nc"

ncform<-c(
   "argo",
   "standard1",
   "standard1",
   "standard1",
   "standard1",
   "standard1",
   "standard1",
   "standard1",
   "standard1",
   "standard1",
   "standard1"
)

countyrs<-c(
   1,
   10,
   10,
   10,
   10,
   10,
   10,
   9,
   5,
   10,
   10
)

modnames<-c(
   "argo",
   "cccma_cgcm3_1",
   "gfdl_cm2_1",
   "giss_model_e_r",
   "iap_fgoals1_0_g",
   "inmcm3_0",
   "miroc3_2_hires",
   "mpi_echam5",
   "ncar_ccsm3_0",
   "ncar_pcm1",
   "ukmo_hadcm3"
)

legnames<-c(
   "ARGO",
   "CCCMA_3_1",
   "GFDL_CM2_1",
   "GISS_ER",
   "FGOALS1_0",
   "INMCM3_0",
   "MIROC3_2",
   "ECHAM5",
   "CCSM3_0",
   "PCM1",
   "HADCM3"
)


#
# Function: getncsample
# This function returns a four dimensional array of netCDF potential temperature data.
# In addition to the nc file and source description, the function also
# receives the bounds for the sample:
#  wb: west bound
#  eb: east bound
#  sb: south bound
#  nb: north bound
#  ub: upper bound (ocean depth - usually zero)
#  lb: lower bound
#  sy: start year
#  cy: count years (number of years to sample)
#

getncsample<-function(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy)
{
 if (wb <= eb){
    ptsamp = selectncsamp(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy)
 } else 
 {
    ptsamp1 = selectncsamp(nc, ncsource, wb, 360, sb, nb, ub, lb, sy, cy)
 ptsamp2 = selectncsamp(nc, ncsource,  0,  eb, sb, nb, ub, lb, sy, cy)
    ptsamp  = abind(ptsamp1, ptsamp2, along=1)
 }
 return(ptsamp)
}

#
# Function: selectncsample
# Called from getncsample, this does the actual data retrieval from
# the nc file.  Same parameters as getncsample
#

selectncsamp <- function(nc, ncsource, wb, eb, sb, nb, ub, lb, sy, cy){

 if (ncsource=="argo"){ 
    longvar<-"LONGITUDE"
    latvar<-"LATITUDE"
    depvar<-"LEVEL"
    timevar<-"TIME"
    tempvar<-"PTEMP"
    k<-0} else
    if (ncsource=="standard1"){ 
       longvar<-"lon"
       latvar<-"lat"
       depvar<-"depth"
       timevar<-"time"
       tempvar<-"thetao"
       k<-273.15}
 
 long<-get.var.ncdf(nc,longvar)
 lat<-get.var.ncdf(nc,latvar)
 depth<-get.var.ncdf(nc,depvar)
 
 longx<-which(long>=wb & long<=eb)
 latx<-which(lat>=sb & lat<=nb)
 depthx<-which(depth>=ub & depth<=lb)
 
 startlong<-min(longx)
 countlong<-length(longx)
 
 startlat<-min(latx)
 countlat<-length(latx)
  
 startdepth<-min(depthx)
 countdepth<-length(depthx)
 
 sm=(sy-1)*12+1
 cm=cy*12
 
 start<-c(startlong,startlat,startdepth,sm)
 count<-c(countlong,countlat,countdepth,cm)
 
 ptsamp<-get.var.ncdf(nc,tempvar,start,count) - k
 dim(ptsamp)<-count

 attributes(ptsamp)$dimnames <- list(long[longx],lat[latx],depth[depthx],1:cm-0.5)

 return(ptsamp)
 
}

#
# Function: sinefit
# This function fits a sine wave to 12 months of data and returns
# the mean, phase, amplitude, and R^2 verification statistic.
#
# The input variable y is a two dimensional array with rows representing depths
# and columns representing months.  The variable phaseoffsetdays is used to mark
# the desired equinox.
#

sinefit<-function(y,phaseoffsetdays)
{
 y<-t(y)
 x<-(2*pi*(((0:11+.5)/12)+(phaseoffsetdays/365)))
 
 x1<-cos(x)
 x2<-sin(x)
 
 m<-lm(y ~ x1 + x2)
 s<-summary(m)
 r<-coef(m)
 a<-r[1,]
 b<- -atan2(r[2,],r[3,])
 b<-ifelse(b<0,2*pi+b,b)
 b<-(b*365)/(2*pi)
 c<-sqrt(r[2,]^2+r[3,]^2)
 wmean<-a
 wphase<-b
 wamp<-c
 wr2<-NULL
 for (i in 1:length(s)) {wr2<-c(wr2,s[[i]]$adj.r.squared)} 
 ret<-cbind(wmean,wphase,wamp,wr2)
 return(ret)
}

#
#### MAIN PROGRAM ####

outlst<-NULL
ptlst<-NULL

date()
  
for (i in 1:length(ncfiles))
{
  #
  # Get the sample of data
  #
  
  ncsource1<-ncform[i]
  cy<-countyrs[i] # number of years

  nc1<-open.ncdf(ncfiles[i])
  
  if (ncsource1=="argo"){ 
     latvar<-"LATITUDE"} else
     if (ncsource1=="standard1"){ 
        latvar<-"lat"}

  nclat<-get.var.ncdf(nc1, latvar)

  lats<-nclat[which(nclat>=sb & nclat<=nb)]

  phl<-NULL
  aml<-NULL
  mnl<-NULL
  r2l<-NULL

  phw<-NULL
  amw<-NULL
  mnw<-NULL
  r2w<-NULL

  wghtv<-NULL
    
  for (j in lats){
  

    #
    # Determine the phase offset days for the desired equinox.
    # Calculations may be off by a day or two.
    # The switch in hemispheres is at the center of the ITCZ ~8N
    #

    if (as.numeric(j)> 7) {pod<- -80.25} else {pod<- 102.25}

    ptnc1<-getncsample(nc1, ncsource1, wb, eb, j, j, ub, lb, sy, cy)  
   
    d<-dim(ptnc1)
    dn<-dimnames(ptnc1)
    dim(ptnc1)<-c(d[1],d[2],d[3],12,cy)
    attributes(ptnc1)$dimnames <- list(dn[1][[1]],dn[2][[1]],dn[3][[1]],0:11+0.5,1:cy)
    ptnc1<-apply(ptnc1,c(1,2,3,4),mean)
  
    #
    # Get rid of longitudes that have NA values.
    #
    longs<-apply(ptnc1,1,function(x)all(!is.na(x)))
    longs<-which(longs==TRUE)
    ptsamp1<-ptnc1[longs,,,,drop=FALSE]

    ptmean1<-apply(ptsamp1[,,,,drop=FALSE],c(3,4),mean)
    ptfit1<-sinefit(ptmean1,pod)
  
    phc<-ptfit1[,"wphase",drop=FALSE]
    amc<-ptfit1[,"wamp",drop=FALSE]
    mnc<-ptfit1[,"wmean",drop=FALSE]
    r2c<-ptfit1[,"wr2",drop=FALSE]
  
    phl<-cbind(phl,phc)
    aml<-cbind(aml,amc)
    mnl<-cbind(mnl,mnc)
    r2l<-cbind(r2l,r2c)

    wght<-length(longs)*cos(2*pi*as.numeric(j)/360)
    wghtv<-c(wghtv,wght)
    
    phw<-cbind(phw,phc*wght)
    amw<-cbind(amw,amc*wght)
    mnw<-cbind(mnw,mnc*wght)
    r2w<-cbind(r2w,r2c*wght)

  }

  dimnames(phl)[[2]]<-lats
  dimnames(aml)[[2]]<-lats
  dimnames(mnl)[[2]]<-lats
  dimnames(r2l)[[2]]<-lats

  phw<-apply(phw,1,sum)/sum(wghtv)
  amw<-apply(amw,1,sum)/sum(wghtv)
  mnw<-apply(mnw,1,sum)/sum(wghtv)
  r2w<-apply(r2w,1,sum)/sum(wghtv)
    
  ptl<-abind(phl,aml,mnl,r2l,along=0,new.names=c("wphase","wamp","wmean","wr2"))
  ptw<-cbind(phw,amw,mnw,r2w)

  dimnames(ptw)[[2]]<-c("wphase","wamp","wmean","wr2")

  tmplst<-list(modname=modnames[i],legname=legnames[i],ncfile=ncfiles[i],ptl=ptl,ptw=ptw)
  ptlst<-c(ptlst,list(tmplst))

  save(ptlst,file=outfile)

  close.ncdf(nc1)

}

date()

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

#
# Model Sine Fit Analysis - Image Plots
#

library("fields")

load("pt-55s-to-55n-mtf.RData")

xout<-seq(0,1000,10)

i<-1

############## AMP ################

mamp <- ptlst[i][[1]]$ptl["wamp",,]
xm <- as.numeric(dimnames(mamp)[[1]])
ym <- as.numeric(dimnames(mamp)[[2]])

spamp<-apply(mamp,2,function(x){spline(xm,x,xout=xout)$y})
pramp<-apply(spamp,1,function(x){ecdf(x)(x)})

dimnames(pramp)[[1]]<-ym
dimnames(pramp)[[2]]<-xout

x<-as.numeric(dimnames(pramp)[[1]])
y<-as.numeric(dimnames(pramp)[[2]])
  
plotname<-"ARGO Amplitude - Percentile Rank by Depth"

dev.new()
image.plot(x, y, pramp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))


png("sine-amp-argo.png", bg="transparent", width=450, height=450)

image.plot(x, y, pramp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

dev.off()

############## R2 #################

mamp <- ptlst[i][[1]]$ptl["wr2",,]
xm <- as.numeric(dimnames(mamp)[[1]])
ym <- as.numeric(dimnames(mamp)[[2]])

spamp<-t(apply(mamp,2,function(x){spline(xm,x,xout=xout)$y}))

dimnames(spamp)[[1]]<-ym
dimnames(spamp)[[2]]<-xout

x<-as.numeric(dimnames(spamp)[[1]])
y<-as.numeric(dimnames(spamp)[[2]])
  
plotname<-"ARGO R2 Stat"

dev.new()
image.plot(x, y, spamp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

png("sine-r2-argo.png", bg="transparent", width=450, height=450)

image.plot(x, y, spamp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

dev.off()

############## PHASE #################

mamp <- t(ptlst[i][[1]]$ptl["wphase",,])

xm <- as.numeric(dimnames(mamp)[[1]])
ym <- as.numeric(dimnames(mamp)[[2]])

plotname<-"ARGO Lag (Days)"
  
dev.new()
image.plot(xm, ym, mamp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

png("sine-phase-argo.png", bg="transparent", width=450, height=450)

image.plot(xm, ym, mamp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

dev.off()

############## MEAN ################

mamp <- ptlst[i][[1]]$ptl["wmean",,]
xm <- as.numeric(dimnames(mamp)[[1]])
ym <- as.numeric(dimnames(mamp)[[2]])

spamp<-apply(mamp,2,function(x){spline(xm,x,xout=xout)$y})
pramp<-apply(spamp,1,function(x){ecdf(x)(x)})

dimnames(pramp)[[1]]<-ym
dimnames(pramp)[[2]]<-xout

x<-as.numeric(dimnames(pramp)[[1]])
y<-as.numeric(dimnames(pramp)[[2]])
  
plotname<-plotname<-"ARGO Mean - Percentile Rank by Depth"
  
dev.new()
image.plot(x, y, pramp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

png("sine-mean-argo.png", bg="transparent", width=450, height=450)

image.plot(x, y, pramp, main=plotname,
              ylab="Depth (M)", xlab="Latitude", ylim=rev(range(y)))

dev.off()

Comments