Bathymetry data/Rcode 2

#grid data in ASCII raster format with no headers from  http://www.ngdc.noaa.gov/mgg/gdas/gd_designagrid.html
  z = read.table('http://wiki.cbr.washington.edu/qerm/images/a/ad/Eli_bath.txt',header=F)
  x = seq(130,180,1/15) #4-minute grid
  y = seq(30,70,1/15)

  z = z[nrow(z):1,]     #why it needs to be flipped over, who knows?
  z <- t(as.matrix(z))  #or transposed for that matter....

# obtain average of four corners of each square ... this is the basis for the color matrix

  z2 <- (z[-1, -1] + z[-1, -ncol(z)] + z[-nrow(z), -1] + z[-nrow(z), -ncol(z)])/4

  zmin=min(z2)
  zmax=max(z2)

  NwaterColors = round(abs(zmin))
  NlandColors = round(abs(zmax))

  #taking first 1/3 of topo.colors of length NwaterColors for water
    water.colors = topo.colors(NwaterColors*3)[1:NwaterColors]
  #taking terrain.colors of length NlandColors for land
    land.colors = terrain.colors(NlandColors)
    all.colors = c(water.colors, land.colors)

    z3 <- z2-min(z2)
    color.matrix <- all.colors[round(z3)] 

# the perspective plot 

  phi <- 60   # height angle
  theta <- 0  # degrees east (0 means looking due north)

  persp(x,y,z,col=color.matrix,
		phi=phi,theta=theta,ltheta=120,shade=0.5,bor=NA,
		 xlab="Longitude", ylab="Latitude", zlab="Elevation (m)",
		expand=0.5)

• This code produces
  
• Return to Bathymetry Data page