sh...@ucar.edu wrote: > I would like to create a 3d topographic map using lat/lon and > z(height). I have been scouring the R help pages and have not located > the package I am looking for. Does anyone have a suggestion of package > that will work for this?
I have some code for generating some pretty cool 3D topographic maps, even interactive ones. It's not very difficult to do, really, and much of the code is actually just for generating nice colours for the map. :-) I thought it was quite surprising and interesting to see how shallow the ocean is many place, and how much the depth actually varies. For instance, the North Sea (between the coast of Norway and the coast of England) is only about 100 meters deep most places. But in other places the ocean is several kilometers deep. Anyway, here's the code: # First download the SRTM30_Plus data file (about 55 MiB) from # ftp://topex.ucsd.edu/pub/srtm30_plus/srtm30/data/w020n90.Bathymetry.srtm # See also http://topex.ucsd.edu/WWW_html/srtm30_plus.html # Bounding box minlat=40 maxlat=90 minlong=-20 maxlong=20 # Range of bounding box longrange=maxlong-minlong latrange=maxlat-minlat # Number of columns and rows of data in the bathymetry file ncols=60 * 2 * longrange nrows=60 * 2 * latrange ntot=ncols*nrows # Total number of points -- probably too many to plot ... steps=10 # ... so we use only every 'steps' point in each direction ux=seq(1,ncols,by=steps) # X coordinates (longitude) uy=seq(1,nrows,by=steps) # Y coordinates (latitude) dat=expand.grid(x=ux,y=uy) # Create grid dat=transform(dat,index=(y-1)*ncols+x) # Calculate byte index in file # Transform indices to real coordinates dat=transform(dat, long=ggplot2::rescale(x, from=c(1,ncols), to=c(minlong+1/60/4, maxlong-1/60/4), clip=FALSE), lat=ggplot2::rescale(y, from=c(1,nrows), to=c(maxlat-1/60/4, minlat+1/60/4), clip=FALSE)) # Name of file containing the bathymetry data filename="w020n90.Bathymetry.srtm" # Read bathymetry data from the file library(R.utils) dat$z=readBinFragments(filename, what="integer", size=2, signed=TRUE, idxs=dat$index,endian="big") # Create matrix containing the data dat.mat=matrix(dat$z,ncol=length(ux),byrow=TRUE) dat.mat=t(dat.mat)[,nrow(dat.mat):1] # Calculate topographic colours for elevation values map.colours=function(z, water=c(-5000,0), land=c(0,1500), mountain,...) { z.min=min(z) z.max=max(z) # Divide the land gradient into two parts land=c(land[1],mean(land),land[2]) # Number of unique water and land values nw=length(unique(z[z<0])) nl=length(unique(z[z>=0])) # Create the output colour component vectors h <- s <- v <- rep(1, length(z)) # h-values for water hvals=c(43/60,31/60) # Create water colours if(nw==1) h[z<0]=hvals[2] else if(nw>1) h[z<0]=approx(water,hvals,xout=z[z<0],rule=2)$y # h, s and v values for land hvals=c(4/12,2/12,0/12) svals=c(1,1,0) vvals=c(.65,.9,.95) # Create land colours if(nl==1) { h[z>=0] = hvals[1]; s[z>=0] = svals[1]; v[z>=0] = vvals[1]} else if(nl>1) { h[z>=0]=approx(land,hvals,xout=z[z>=0],rule=2)$y s[z>=0]=approx(land,svals,xout=z[z>=0],rule=2)$y v[z>=0]=approx(land,vvals,xout=z[z>=0],rule=2)$y } hsv(h,s,v) } # Basic contour plot -- not very pretty ... with(dat, contour(unique(long),rev(unique(lat)),dat.mat)) # Colours make the whole thing look quite pretty ... ncolours=100 image(unique(dat$long),rev(unique(dat$lat)), matrix(1:length(dat.mat),nrow=nrow(dat.mat)), col=map.colours(dat.mat),xlab="Longitude",ylab="Latitude") # And it looks even prettier in 3D ... zfacet <- dat.mat[-nrow(dat.mat), -ncol(dat.mat)] with(dat, persp(unique(long), rev(unique(lat)), dat.mat, expand=.2, col=map.colours(zfacet), border=NA, shade=.3, phi=60, theta=0, axes=FALSE)) # How about an *interactive* 3D map ... # Try holding the left mouse button and dragging # Now try the mouse wheel (or right button and dragging) # Finally, try clicking the middle button (or wheel) and dragging library(rgl) with(dat, surface3d(unique(long), rev(unique(lat)), dat.mat/1000, col=map.colours(dat.mat))) # Note that the elevation values are of course very much exaggerated # In 'real life', it would look something like this (i.e., completely flat) with(dat, surface3d(unique(long), rev(unique(lat)), dat.mat/1000000, col=map.colours(dat.mat))) -- Karl Ove Hufthammer ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.