[R] Oceanographic lattice plots

2003-10-19 Thread Sam McClatchie
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 - 

Re: [R] Oceanographic lattice plots?

2003-10-18 Thread Deepayan Sarkar
On Saturday 18 October 2003 13:14, Mike Prager wrote:
 R 1.8.0 on Windows XP Professional.  A huge THANK YOU to the R Team for
 this marvelous software.

 I am making lattice plots of oceanographic data.  The usual layout does not
 conform to plotting conventions that marine scientists use when depth is
 the independent variable.  Under those conventions, plots are made with the
 origin at the upper left, depth on the vertical axis (increasing as it goes
 down), and the dependent variable on the horizontal axis (increasing to the
 right).

In case you decide to work on this yourself, this might be useful: 

Ideally, specifications like 

xyplot(depth ~ x, ylim = c(10, 0)) 

should reverse the y-axis direction. As pointed out some time back, this 
doesn't work in lattice currently, but that should be fixed in the future (it 
already works in my development version).

 That convention has implications not just in how axes are labeled and set
 up, but also when using smoothing routines such as panel.lowess(), because
 the smoothed values are on the horizontal axis, not the vertical axis.

These should be easy to modify.

Deepayan

 Before I start looking at and modifying the R code that makes up the
 relevant routines, I wonder if any reader has already developed R routines
 for this purpose?

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help