Bathymetry data/Rcode 2

From QERM Wiki
Jump to: navigation, search

#grid data in ASCII raster format with no headers from
  z = read.table('',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


  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)

		 xlab="Longitude", ylab="Latitude", zlab="Elevation (m)",
Personal tools