makeTopography: Make Topography Data From Freely Available Online Data

Description

Make topography data suitable for the contour and contourLines functions using freely available global seafloor topography data.

Usage

makeTopography (dat, digits=2, func=NULL)

Arguments

dat

data frame with three optionally-named columns: X, Y, and Z. The columns must appear in that order.

digits

integer indicating the precision to be used by the function round on (X,Y) values.

func

function to summarize Z if (X,Y) points are duplicated. Defaults to mean() if no function is specified.

Value

List with elements x, y, and z. x and y are vectors, while z is a matrix with rownames x and colnames y. contour and contourLines expect data conforming to this list format.

Details

Data obtained through the acquisition form at http://topex.ucsd.edu/cgi-bin/get_data.cgi is suitable for this function. read.table will import its ASCII files into R/S, creating the data argument for this function.

When creating data for regions with longitude values spanning -180\(^\circ\) to 0\(^\circ\), consider subtracting 360 from the result's X coordinates (x).

When creating bathymetry data, consider negating the result's elevations (z) to give depths positive values.

Combinations of (X,Y) do not need to be complete (z[x,y]=NA) or unique (z[x,y]=func(Z[x,y])).

Examples

Run this code

local(envir=.PBSmapEnv,expr={
  oldpar = par(no.readonly=TRUE)
  #--- Example 1: Sample data frame and conversion.
  file <- data.frame(X=c(1,1,2,2),Y=c(3,4,3,4),Z=c(5,6,7,8))
  print(makeTopography(file))

  #--- Example 2: Aleutian Islands bathymetry
  isob <- c(100,500,1000,2500,5000)
  icol <- rgb(0,0,seq(255,100,len=length(isob)),max=255)
  afile <- paste(system.file(package="PBSmapping"),
    "/Extra/aleutian.txt",sep="")
  aleutian <- read.table(afile, header=FALSE, col.names=c("x","y","z"))
  aleutian$x <- aleutian$x - 360
  aleutian$z <- -aleutian$z
  alBathy <- makeTopography(aleutian)
  alCL <- contourLines(alBathy,levels=isob)
  alCP <- convCP(alCL)
  alPoly <- alCP$PolySet
  attr(alPoly,"projection") <- "LL"
  plotMap(alPoly,type="n")
  addLines(alPoly,col=icol)
  data(nepacLL,envir=.PBSmapEnv)
  addPolys(nepacLL,col="gold")
  legend(x="topleft",bty="n",col=icol,lwd=2,legend=as.character(isob))
  par(oldpar)
})