[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