Bathymetry data/Rcode 2
From QERM Wiki
(Difference between revisions)
m |
m (fixed category) |
||
Line 42: | Line 42: | ||
* '''[[Bathymetry Data | Return to Bathymetry Data page]]''' | * '''[[Bathymetry Data | Return to Bathymetry Data page]]''' | ||
− | [[Category: R | + | [[Category: R tips]] |
Latest revision as of 23:33, 21 July 2010
#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