[R] Oceanographic lattice plots
Sam McClatchie
s.mcclatchie at niwa.co.nz
Sun Oct 19 23:13:44 CEST 2003
System info:
Red Hat 9.0
R Version 1.8.0
ESS 5.1.21
Emacs 21.2.1
-------------------
Hello
I've been working with Paul Murrell here in New Zealand to develop plots
of temperature and density profiles at 22 stations spanning a frontal
zxone, and then overlaying the isothermal and mixed layer depths
(similar to Kara et al. 2003 JGR Oceans 108(C3) pg. 24-5, figure 2).
Paul Murrell has a package called gridBase which works with R 1.8.0 that
does the job nicely.
I've put the plot on a publicly accessible ftp site
<ftp.niwa.co.nz/incoming/mcclatchie/misc/ocean.plot.pdf>
Here is the code used to generate the plot that Paul developed with my
collaboration.
---------------------------------------
panel.t.s.profiles.mld.paul.alternate <- function()
{
## Purpose:
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Paul Murrell & Sam McClatchie, Date: 16 Oct 2003, 09:43
require(gridBase)
## zt, p, zsigma, and lat are in the workspace
## ILD data generated by "calculate.MLD.and.ILD(0.8)" etc
#ILD values are in the workspace
#convert to positive values
ILD0.2 <- -ILD0.2
ILD0.2 <- matrix(c(ILD0.2,NA,NA),ncol=8,byrow=T)
ILD0.2.vec <- c(ILD0.2[1,],ILD0.2[2,],ILD0.2[2,])
ILD0.5 <- -ILD0.5
ILD0.5 <- matrix(c(ILD0.5,NA,NA),ncol=8,byrow=T)
ILD0.5.vec <- c(ILD0.5[1,],ILD0.5[2,],ILD0.5[2,])
ILD0.8 <- -ILD0.8
ILD0.8 <- matrix(c(ILD0.8,NA,NA),ncol=8,byrow=T)
ILD0.8.vec <- c(ILD0.8[1,],ILD0.8[2,],ILD0.8[2,])
## dummy MLD data
MLD0.5 <- ILD0.5*0.75
MLD0.5.vec <- c(MLD0.5[1,],MLD0.5[2,],MLD0.5[2,])
nrow <- 3
ncol <- 8
grid.newpage()
par(xpd=NA, col.axis="grey")
# Use grid to make a layout of rows and columns
# Each "row" consists of a plot region with a label area on top
# (i.e., a Trellis-like arrangement) hence the "nrow*2".
# The label area is 1 line high, the plot areas
# consume the remaining height.
# Maybe you could add extra rows and cols to this layout to
# create small gaps between each plot
# This layout is within a viewport which leaves margins for axes
# and labels
push.viewport(viewport(x=unit(4, "lines"),
y=unit(4, "lines"),
width=unit(1, "npc") - unit(6, "lines"),
height=unit(1, "npc") - unit(6, "lines"),
just=c("left", "bottom"),
layout=grid.layout(nrow*2, ncol,
heights=unit(rep(c(1, 1), nrow),
rep(c("lines", "null"), nrow)))))
for (i in 1:nrow) {
for (j in 1:ncol) {
index <- (i-1)*ncol + j
evenrow <- i %% 2 == 0
evencol <- j %% 2 == 0
if (index < 23) {
# Go to plot region i, j
# Set the yscale for doing the overlay later
push.viewport(viewport(layout.pos.row=i*2, layout.pos.col=j,
yscale=c(400, 0)))
grid.rect(gp=gpar(col="grey"))
# Draw first plot
# Here's where we use gridBase to put a plot into a grid viewport
# The par(plt=gridPLT()) makes the plotting region line up with
# the current grid viewport (pushed two lines ago)
par(plt=gridPLT(), new=TRUE, cex=0.8)
plot(zt[,index],-p[,1], ylim=c(400,0), xlim=c(6,15), type='l',
xlab="",ylab="", axes=FALSE, xpd=FALSE)
# Draw axes (only do some)
if (j == 1 && !evenrow)
axis(2, col="grey")
if (i == nrow && !evencol)
axis(1, col="grey")
par(new=TRUE)
# Draw second plot
plot(zsigmat[,index],-p[,1], ylim=c(400,0), type='l',lty=2,
xlab="",
ylab="",
xlim=c(26.1,27.4),
axes=FALSE, xpd=FALSE)
# Draw axes
if ((j == ncol || index == 22) && evenrow)
axis(4, col="grey")
pop.viewport()
# Draw the latitude labels
push.viewport(viewport(layout.pos.row=i*2 - 1, layout.pos.col=j,
gp=gpar(col="grey", fill="light grey")))
grid.rect()
grid.text(round(lat[index,],digits=2), gp=gpar(col="white"))
# Draw top axes
if (i == 1 && evencol) {
par(plt=gridPLT())
axis(3, col="grey")
}
pop.viewport()
}
}
}
# Overlay mixed layer depths lines
# Here we use some grid functions to do drawing
# The 0.5 means "half way across the region",
# the "native" means that the value is relative
# to the yscale we set up when we created the viewport
for (i in 1:nrow) {
for (j in 1:ncol) {
index <- (i-1)*ncol + j
if (index < 23) {
push.viewport(viewport(layout.pos.row=i*2, layout.pos.col=j,
yscale=c(400, 0)))
# Do a move.to in the first column and a line.to otherwise
if (j == 1)
grid.move.to(0.5, unit(MLD0.5[i, j], "native"))
else
grid.line.to(0.5, unit(MLD0.5[i, j], "native"),
gp=gpar(lty="dotted"))
pop.viewport()
}
}
}
# Draw the overlay points last to overwrite the overlay lines
for (i in 1:nrow) {
for (j in 1:ncol) {
index <- (i-1)*ncol + j
if (index < 23) {
push.viewport(viewport(layout.pos.row=i*2, layout.pos.col=j,
yscale=c(400, 0)))
# Overlay mixed layer depths points
grid.points(0.5, unit(ILD0.5[i, j], "native"), pch=21,
size=unit(3, "mm"),
gp=gpar(fill=NULL))
grid.points(0.5, unit(ILD0.2[i, j], "native"), pch=21,
size=unit(3, "mm"),
gp=gpar(fill=NULL))
grid.points(0.5, unit(ILD0.8[i, j], "native"), pch=21,
size=unit(3, "mm"),
gp=gpar(fill=NULL))
grid.points(0.5, unit(MLD0.5[i, j], "native"), pch=21,
size=unit(3, "mm"),
gp=gpar(fill="grey"))
pop.viewport()
}
}
}
# Do some overall axis labels
push.viewport(viewport(layout.pos.row=nrow*2))
grid.text(expression(paste("temp", .^o,"C", sep="")),
y=unit(-3, "lines"))
pop.viewport()
push.viewport(viewport(layout.pos.col=1))
grid.text("depth m", x=unit(-3, "lines"), rot=90)
pop.viewport()
pop.viewport()
}
------------------
I hope that some of this may be of use to you. It would be nice if you
could let both myself and Paul see the result of efforts.
Cheers
Sam
--
Sam McClatchie, Research scientist (fisheries acoustics))))))))))
NIWA (National Institute of Water & Atmospheric Research Ltd)
PO Box 14 901, Kilbirnie, Wellington, New Zealand
s.mcclatchie at niwa.cri.nz
Research home page <http://www.smcc.150m.com/>
/\
>><xX(&>
/// \\\
//// \\\\
/// <%)Xx><<
///// \\\\\\
><(((@>
><(((%> ..>><xX(?>O<?)Xx><<
More information about the R-help
mailing list