## Computer code and used to generate the results for the Journal of Climate paper "An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity" by Nicholas Lewis, Bath, uk




Download 0.75 Mb.
Name## Computer code and used to generate the results for the Journal of Climate paper "An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity" by Nicholas Lewis, Bath, uk
page1/13
A typeDocumentation
manual-guide.com > manual > Documentation
  1   2   3   4   5   6   7   8   9   ...   13
## Computer code and used to generate the results for the Journal of Climate paper "An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity" by Nicholas Lewis, Bath, UK

## objectivef06-code-jclim.doc Copyright Nicholas Lewis May 2013

### R code for working on the Forest 2006 data and generating climate sensitivity PDFs etc, and example results

# The code has been run and tested on a Windows 7 64-bit system with 8Gb memory.

# Acronyms used: F06 = Forest et al. 2006; CSF05 or Curry = Curry, Sanso and Forest 2005; SFZ08 or Sanso = Sanso, Forest and Zantedeschi 2008

#####################################################

### Set up the R workspace and load the required data

#####################################################

### Open a clean R workspace.

### Manually install (from a CRAN mirror site) and load the following required R packages: ncdf, lattice, Rcolorbrewer, akima and fields

### Create & set working directory to desired path (sub-directories 'data', in which all data is stored, 'plots' in which plots are saved and 'code', in which code can be saved, are also created). A different path can be used.

# Example path - alter as necessary to suit the system disk and directory structure

path="C:/R/Forest06/fin"

dir.create(path)

setwd(path)

getwd()

dir.create(paste(path, "/data", sep=""))

dir.create(paste(path, "/data/orig", sep=""))

dir.create(paste(path, "/plots", sep=""))

dir.create(paste(path, "/code", sep=""))

dir.create(paste(path, "/GRL06_reproduce", sep=""))

dir.create(paste(path, "/GRL06_reproduce/data", sep=""))

##########################################################################################################

### Download the following archive files and extract the stated data files to the indicated sub-directory:

##########################################################################################################

## Download the archive file 'brunodata_May23.zip'

# extract the only file the zip archive contains, brunodata_May23.tgz, to the path/data/orig directory

# brunodata_May23.tgz contains data used in Sanso, Forest and Zantedeschi (2008) 'Inferring Climate System Properties Using a Computer Model' (SFZ08), supplied by Chris Forest to Bruno Sanso (kindly provided by Sanso). This data has been confirmed by Chris Forest to be the same as that actually used in F06. There are two sets of surface data, identical apart from the surface diagnostic control data matrix, the one derived from the HadCM2 model (subsequently extracted and stored as sfc.ctl.sH) being used, as if F06. There was no deep ocean data in the brunodata_May23.tgz archive.

## download the GRL06_reproduce4.zip archive file

# This file contain five files made into a zip file after extracting them from the data directory of GRL06_reproduce.tgz, a 2 GB archive file at http://svante.mit.edu/research/IGSM/data/IGSM_1/GRL06_reproduce.tgz

# extract the Levitus 2005 pentadal 0-3000m heat content data file, the HadCRUT observational and the HadCM2 control run decadal mean surface temperature data file and the ct2.vl4 upper air raw control data file from the GRL06_reproduce4.zip archive into the path/GRL06_reproduce/data directory

# hc5yr-w0-3000m.dat

# obs_dec_9_1896-9_1996.nc

# hadcm2.nc

# gfdl.ctrl.r30.1000.dat

# ct2.vl4

## Download the Lewis2013data.zip archive file

#extract the following files from Lewis2013data.zip to the path/data/orig directory:

# HadCM3_vol_mean_0-3039m.txt

# hadcrut4_median.nc

# F06_08.pdf.S.Rd

#These files contain HadCM3 6100 year control run deep ocean temperature data (HadCM3_vol_mean_0-3039m.txt); HadCRUT4 surface gridded monthly temperature observational data (hadcrut4_median.nc, downloaded 16 April 2012); digitised data from main results climate sensitivity PDFs per figures in F06 and F08 (F06_08.pdf.S.Rd).

# extract the following files from the Lewis2013data.zip into the path/GRL06_reproduce/data directory.

# dto.gsovsv.datGSOVSV

# mit.mdl.txt

#the dto.gsovsv.datGSOVSV file is MIT model 0-3000m deep ocean temperature data The file is produced by the following line in pro dto when running the GRL06_reproduce IDL code: # openw,1,'~/GRL06_reproduce/processed_gsovsv/dto.gsovsv.dat'+runlist

# the mit.mdl.txt file is generated from running a modified version of the IDL code in the GRL06_reproduce.tgz archive that generates a text file containing y/e 31Aug1861 to y/e 31Aug2001 annual mean 4 degree latitude band surface temperature MIT model data.

# The very large HadCM2 surface and upper air annual mean control data (sfc999, sfc1, sfc691 and ua.900) is in a separate archive file HadCM2.zip and is not essential.

# If wished, download the HadCM2.zip file and extract the following files to the path/data/orig directory:

# sfc999

# sfc1

# sfc691

# ua.900

## Ignore the lines, which are for information only

#Portion of mit2d_diag.pro code:

#; corrected to GSOVSV time period 1995 = yr 134

#t_mit = findgen(135)+1861. ; 1-SEP-1860 to 31-AUG-1995 - per C E Forest email model y/e was 30 Nov not 31 Aug

#Sample journal output when IDL code is running:

#;MIT run, rrr, yy:101.03 101 03

#;running ~/GRL06_reproduce/antsrfzonall

#;anexp, year1, year2, nyt, dty, jm 101.030

#; 1861 2001 141

#; 1.00000

#; 46

########################################################

### Code functions used to perform the main computations

########################################################

### Input all the below code functions, by copying and pasting all the code to the R console

### Utility functions - purpose as described in first line of each

ssq = function(x, na.rm=FALSE) sum(x*x, na.rm=na.rm)# computes sum of squares, ignoring NA values if na.rm=TRUE

repmat= function(X,m,n=1){ # tiles matrix X, m rowwise by n columnwise times, as in Matlab, treating a vector as a 1 row (not 1 column) matrix

if( is.vector(X) ) { X= matrix(X, nrow=1) }

mx = dim(X)[1]

nx = dim(X)[2]

matrix(t(matrix(X,mx,nx*n)),mx*m,nx*n,byrow=TRUE)

}

resetPar= function() { # resets graphical parameters to default values
    dev.new()
    op= par(no.readonly = TRUE)
    dev.off()
    op
}

position= function(x, dims=c(101,81,41)) { # expresses a linear position x as the dimensioned position in an array of dimension dims

test= array(0, dim=dims)

test[x]= 1

posn=vector()

for ( i in 1:length(dims) ) { posn[i]= which( apply(test, i, sum)==1 ) }

posn

}

positions= function(x, dims=c(101,81,41)) { # expresses a vector x of linear positions as the dimensioned positions in an array of dimension dims

n= length(x)

posn=matrix(nrow=n, ncol=3)

posn[,3]= 1 + (x-1) %/% (dims[1] * dims[2])

x= (x-1) %% (dims[1] * dims[2])

posn[,2]= 1 + x %/% dims[1]

posn[,1]= 1 + x %% dims[1]

posn

}

vecPos= function(x, y, z, dims=c(101,81,41)) { # expresses a 3D coordinate position x, y, z in an array of dimension dims as a vectorised position

x + (y-1) * dims[1] + (z-1) * dims[1] * dims[2]

}

deltaCoords= function(coords) { # gives coordinates of adjacent points to a given point's coordinates with one at a time changed, first +, then -

m= length(coords)

deltaCoords= matrix(coords, nrow=2*m, ncol=m, byrow=TRUE)

delta= diag(m)

delta= rbind(delta, -delta)

deltaCoords + delta

}

save(ssq, repmat, resetPar, position, positions, vecPos, deltaCoords, file='code/utilities.Rd')

# Function to compute 5%, 25%, 50%, 75% and 95% CDF points for a matrix whose columns are PDFs, and to add box plots showing these CDF points to the current graph

plotBoxCIs= function(pdfsToPlot, divs, lower, upper, boxsize=0.75, col, yOffset=0.05, spacing=0.75, lwd=3, whisklty=1, plot=TRUE, points=c(0.05,0.25,0.5,0.75,0.95)) {

range= upper - lower

cases= ncol(pdfsToPlot)

boxsize= boxsize

box=boxplot.stats(runif(100,0,1))

stats=matrix(NA,cases,5);

if( identical(points, c(0.05,0.25,0.5,0.75,0.95)) ) { colnames(stats)= c('5%', '25%',' 50%', '75%', '95%') } else { colnames(stats)= as.character(points) }

for (j in 1:cases) {

z= cumsum(pdfsToPlot[,j])

stats[j,1]= which(z > points[1]/divs)[1]*divs + lower - divs

stats[j,2]= which(z > points[2]/divs)[1]*divs + lower - divs

stats[j,3]= which(z > points[3]/divs)[1]*divs + lower - divs

stats[j,4]= which(z > points[4]/divs)[1]*divs + lower - divs

stats[j,5]= which(z > points[5]/divs)[1]*divs + lower - divs

box$stats=matrix(stats[j,], nrow=5)

if( plot==TRUE) { bxp(box, width = NULL, varwidth = FALSE, notch = FALSE, outline = FALSE, names="", plot = TRUE, border = col[j], col = NULL, pars = list(xaxs='i', boxlty=1, boxwex=boxsize, lwd=lwd, staplewex=0.5, outwex=0.5, whisklty=whisklty[j], medlwd=lwd), horizontal= TRUE, add= TRUE, at=boxsize*yOffset+j*boxsize*spacing, axes=FALSE) }

}

stats

}

# function replacing the internal R legend function, which does not implement control over the title font size

legend= function (x, y = NULL, legend, fill = NULL, col = par("col"),

border = "black", lty, lwd, pch, angle = 45, density = NULL,

bty = "o", bg = par("bg"), box.lwd = par("lwd"), box.lty = par("lty"),

box.col = par("fg"), pt.bg = NA, cex = 1, pt.cex = cex, pt.lwd = lwd,

xjust = 0, yjust = 1, x.intersp = 1, y.intersp = 1, adj = c(0,

0.5), text.width = NULL, text.col = par("col"), merge = do.lines &&

has.pch, trace = FALSE, plot = TRUE, ncol = 1, horiz = FALSE,

title = NULL, inset = 0, xpd, title.col = text.col, title.adj = 0.5, title.cex=cex,

seg.len = 2)

{

if (missing(legend) && !missing(y) && (is.character(y) ||

is.expression(y))) {

legend <- y

y <- NULL

}

mfill <- !missing(fill) || !missing(density)

if (!missing(xpd)) {

op <- par("xpd")

on.exit(par(xpd = op))

par(xpd = xpd)

}

title <- as.graphicsAnnot(title)

if (length(title) > 1)

stop("invalid title")

legend <- as.graphicsAnnot(legend)

n.leg <- if (is.call(legend))

1

else length(legend)

if (n.leg == 0)

stop("'legend' is of length 0")

auto <- if (is.character(x))

match.arg(x, c("bottomright", "bottom", "bottomleft",

"left", "topleft", "top", "topright", "right", "center"))

else NA

if (is.na(auto)) {

xy <- xy.coords(x, y)

x <- xy$x

y <- xy$y

nx <- length(x)

if (nx < 1 || nx > 2)

stop("invalid coordinate lengths")

}

else nx <- 0

xlog <- par("xlog")

ylog <- par("ylog")

rect2 <- function(left, top, dx, dy, density = NULL, angle,

...) {

r <- left + dx

if (xlog) {

left <- 10^left

r <- 10^r

}

b <- top - dy

if (ylog) {

top <- 10^top

b <- 10^b

}

rect(left, top, r, b, angle = angle, density = density,

...)

}

segments2 <- function(x1, y1, dx, dy, ...) {

x2 <- x1 + dx

if (xlog) {

x1 <- 10^x1

x2 <- 10^x2

}

y2 <- y1 + dy

if (ylog) {

y1 <- 10^y1

y2 <- 10^y2

}

segments(x1, y1, x2, y2, ...)

}

points2 <- function(x, y, ...) {

if (xlog)

x <- 10^x

if (ylog)

y <- 10^y

points(x, y, ...)

}

text2 <- function(x, y, ...) {

if (xlog)

x <- 10^x

if (ylog)

y <- 10^y

text(x, y, ...)

}

if (trace)

catn <- function(...) do.call("cat", c(lapply(list(...),

formatC), list("\n")))

cin <- par("cin")

Cex <- cex * par("cex")

if (is.null(text.width))

text.width <- max(abs(strwidth(legend, units = "user",

cex = cex)))

else if (!is.numeric(text.width) || text.width < 0)

stop("'text.width' must be numeric, >= 0")

xc <- Cex * xinch(cin[1L], warn.log = FALSE)

yc <- Cex * yinch(cin[2L], warn.log = FALSE)

if (xc < 0)

text.width <- -text.width

xchar <- xc

xextra <- 0

yextra <- yc * (y.intersp - 1)

ymax <- yc * max(1, strheight(legend, units = "user", cex = cex)/yc)

ychar <- yextra + ymax

if (trace)

catn(" xchar=", xchar, "; (yextra,ychar)=", c(yextra,

ychar))

if (mfill) {

xbox <- xc * 0.8

ybox <- yc * 0.5

dx.fill <- xbox

}

do.lines <- (!missing(lty) && (is.character(lty) || any(lty >

0))) || !missing(lwd)

n.legpercol <- if (horiz) {

if (ncol != 1)

warning("horizontal specification overrides: Number of columns := ",

n.leg)

ncol <- n.leg

1

}

else ceiling(n.leg/ncol)

has.pch <- !missing(pch) && length(pch) > 0

if (do.lines) {

x.off <- if (merge)

-0.7

else 0

}

else if (merge)

warning("'merge = TRUE' has no effect when no line segments are drawn")

if (has.pch) {

if (is.character(pch) && !is.na(pch[1L]) && nchar(pch[1L],

type = "c") > 1) {

if (length(pch) > 1)

warning("not using pch[2..] since pch[1L] has multiple chars")

np <- nchar(pch[1L], type = "c")

pch <- substr(rep.int(pch[1L], np), 1L:np, 1L:np)

}

}

if (is.na(auto)) {

if (xlog)

x <- log10(x)

if (ylog)

y <- log10(y)

}

if (nx == 2) {

x <- sort(x)

y <- sort(y)

left <- x[1L]

top <- y[2L]

w <- diff(x)

h <- diff(y)

w0 <- w/ncol

x <- mean(x)

y <- mean(y)

if (missing(xjust))

xjust <- 0.5

if (missing(yjust))

yjust <- 0.5

}

else {

h <- (n.legpercol + (!is.null(title))) * ychar + yc

w0 <- text.width + (x.intersp + 1) * xchar

if (mfill)

w0 <- w0 + dx.fill

if (do.lines)

w0 <- w0 + (seg.len + x.off) * xchar

w <- ncol * w0 + 0.5 * xchar

if (!is.null(title) && (abs(tw <- strwidth(title, units = "user",

cex = cex) + 0.5 * xchar)) > abs(w)) {

xextra <- (tw - w)/2

w <- tw

}

if (is.na(auto)) {

left <- x - xjust * w

top <- y + (1 - yjust) * h

}

else {

usr <- par("usr")

inset <- rep(inset, length.out = 2)

insetx <- inset[1L] * (usr[2L] - usr[1L])

left <- switch(auto, bottomright = , topright = ,

right = usr[2L] - w - insetx, bottomleft = ,

left = , topleft = usr[1L] + insetx, bottom = ,

top = , center = (usr[1L] + usr[2L] - w)/2)

insety <- inset[2L] * (usr[4L] - usr[3L])

top <- switch(auto, bottomright = , bottom = , bottomleft = usr[3L] +

h + insety, topleft = , top = , topright = usr[4L] -

insety, left = , right = , center = (usr[3L] +

usr[4L] + h)/2)

}

}

if (plot && bty != "n") {

if (trace)

catn(" rect2(", left, ",", top, ", w=", w, ", h=",

h, ", ...)", sep = "")

rect2(left, top, dx = w, dy = h, col = bg, density = NULL,

lwd = box.lwd, lty = box.lty, border = box.col)

}

xt <- left + xchar + xextra + (w0 * rep.int(0:(ncol - 1),

rep.int(n.legpercol, ncol)))[1L:n.leg]

yt <- top - 0.5 * yextra - ymax - (rep.int(1L:n.legpercol,

ncol)[1L:n.leg] - 1 + (!is.null(title))) * ychar

if (mfill) {

if (plot) {

fill <- rep(fill, length.out = n.leg)

rect2(left = xt, top = yt + ybox/2, dx = xbox, dy = ybox,

col = fill, density = density, angle = angle,

border = border)

}

xt <- xt + dx.fill

}

if (plot && (has.pch || do.lines))

col <- rep(col, length.out = n.leg)

if (missing(lwd))

lwd <- par("lwd")

if (do.lines) {

if (missing(lty))

lty <- 1

lty <- rep(lty, length.out = n.leg)

lwd <- rep(lwd, length.out = n.leg)

ok.l <- !is.na(lty) & (is.character(lty) | lty > 0)

if (trace)

catn(" segments2(", xt[ok.l] + x.off * xchar, ",",

yt[ok.l], ", dx=", seg.len * xchar, ", dy=0, ...)")

if (plot)

segments2(xt[ok.l] + x.off * xchar, yt[ok.l], dx = seg.len *

xchar, dy = 0, lty = lty[ok.l], lwd = lwd[ok.l],

col = col[ok.l])

xt <- xt + (seg.len + x.off) * xchar

}

if (has.pch) {

pch <- rep(pch, length.out = n.leg)

pt.bg <- rep(pt.bg, length.out = n.leg)

pt.cex <- rep(pt.cex, length.out = n.leg)

pt.lwd <- rep(pt.lwd, length.out = n.leg)

ok <- !is.na(pch) & (is.character(pch) | pch >= 0)

x1 <- (if (merge && do.lines)

xt - (seg.len/2) * xchar

else xt)[ok]

y1 <- yt[ok]

if (trace)

catn(" points2(", x1, ",", y1, ", pch=", pch[ok],

", ...)")

if (plot)

points2(x1, y1, pch = pch[ok], col = col[ok], cex = pt.cex[ok],

bg = pt.bg[ok], lwd = pt.lwd[ok])

}

xt <- xt + x.intersp * xchar

if (plot) {

if (!is.null(title))

text2(left + w * title.adj, top - ymax, labels = title,

adj = c(title.adj, 0), cex = title.cex, col = title.col)

text2(xt, yt, labels = legend, adj = adj, cex = cex,

col = text.col)

}

invisible(list(rect = list(w = w, h = h, left = left, top = top),

text = list(x = xt, y = yt)))

}
  1   2   3   4   5   6   7   8   9   ...   13

Share in:

Related:

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconSession 1 Introducing climate change at as

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconActive climate control systems

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconClimate/terrain: Any Arid / Desert

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconActive climate control systems

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconWeb Extras for Chapter 12: Arid Climate Landscapes

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconThe victorian climate action calendar 5 April – 14 June 2015

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconClimate Classification is an Important Factor in Assessing Quality-of-Care Across Hospitals

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconA guide Specification for the Engineer on the Uponor Climate Cŏntrol™ Network System

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconClimate chnger hot Gas reheat split system air handler

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconBayer confirms leading international position in climate protection...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconSupporting Information for ‘Postponing emission reductions from 2020...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconVirtual compliance to Kyoto by non-state actors and through domestic...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconAbstract Institutional barriers and bridges to local climate change...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconThe following is a step by step detailed instruction set for changing...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconThis document contains a list of references to publications and reports...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconAbstract : Floods induced by either natural events or technical failures...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconHeidelberg launches Prinect Inpress Control 2 – with faster setup...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconGeneral Chipset Sirf starⅢ Frequency L1, 1575. 42 Mhz C/A code 023...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconAbstract: As a forensic consultant and dna fingerprint expert, I...

## Computer code and used to generate the results for the Journal of Climate paper \"An objective Bayesian, improved approch for applying optimal fingerprint techniques to estimate climate sensitivity\" by Nicholas Lewis, Bath, uk iconAbstract fir filter is an important part of the digital filter. The...




manual




When copying material provide a link © 2017
contacts
manual-guide.com
search