Re: [R] 3d topographic map
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[z0]))
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[z0]=hvals[2] else if(nw1)
h[z0]=approx(water,hvals,xout=z[z0],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(nl1)
{
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/100,
col=map.colours(dat.mat)))
--
Karl Ove Hufthammer
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE
Re: [R] 3d topographic map
This is an example that I have found to be very useful example, and one that I have adapted myself: http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=1 On Sun, Jul 25, 2010 at 6:27 PM, David Winsemius dwinsem...@comcast.netwrote: On Jul 25, 2010, at 6:31 PM, sh...@ucar.edu wrote: Hi All- 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? thanks- I suspect most viewers of this message are going to be puzzled. If wireframe and contourplot are not doing it for you, what are the problems? Surely you found references to Lattice and wireframe. The Lattice system has a worked example in Figure 6.11 from this page: http://lmdvr.r-forge.r-project.org/figures/figures.html A search for topographic in one of the r searching facilities is sure to bring multiple hits, . including the usual starting point: https://svn.r-project.org/R/trunk/src/library/graphics/demo/graphics.R So there must be details and issues that you have not chosen to disclose. -- David. __ 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. [[alternative HTML version deleted]] __ 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.
Re: [R] 3d topographic map [SEC=UNCLASSIFIED]
Hi Sherri,
There are examples of topographic maps which you have been pointed to, however,
I suspect that you want to know where you can obtain topographic data from
rather than a canned example.
There are quite a few intricacies to the process so I will go through them for
you.
(1)
Topography files can be found in the geomapdata library. You will probably want
to use the maps package too (if you want to display coastlines). If you are
zoomed in on a certain area, you will probably want to use the hi res coastline
data found in the mapdata package. In the following code, I have also used the
mapproj library which provides the map.grid function for drawing lat/lon lines
on map projections.
rm(list=ls())
library(maps)
library(mapdata)
library(geomapdata)
library(mapproj)
(2) Load the variables
data(ETOPO5)
xf - attr(ETOPO5,'lon') # all lon
yf - attr(ETOPO5,'lat') # all lat
zf - ETOPO5 # elevation
(3) The topography data has to be flipped horizontally and vertically and you
also require the data to be placed on a regular grid and you need to specify
limits,
xn - which(xf=your_lon_low_limit xf=your_lon_high_limit)
yn - which(yf=your_lat_low_limit yf=your_lat_high_limit )
x - xf[xn]
y - yf[yn]
z - zf[xn,yn]
za - apply(z,1,rev)
zb - t(za)
yb - -rev(y)
xb - x
nx - length(xb)
dx - (xb[nx]-xb[1])/(nx-1)
ny - length(yb)
dy - (yb[ny]-yb[1])/(ny-1)
mf - 5
xo - seq(xb[1],xb[nx],by=dx/mf)
yo - seq(yb[1],yb[ny],by=dy/mf)
nxo - length(xo)
nyo - length(yo)
zo1 - matrix(NA,nrow=nxo,ncol=nyo)
(4) Fill your matrix with the data
io - 1
for (i in 1:nx) {
jo - 1
for (j in 1:ny) {
zo1[io,jo] - zb[i,j]
jo - jo+mf
}
io - io+mf
}
(5) Smooth the image
zo - image.smooth(zo1,theta=1.75) (you should try different values of theta)
image(xo,yo,zo$z) will produce a topographic map for you.
contour(xo,yo,zo$z) will produce contour plots of the data for you (which you
can overlay if you like).
map('worldHires',...) will draw coastlines for you
Hope this helps.
Cheers,
Justin
Justin Peter
Research Scientist
Earth System Modelling and Radar Applications Group
Centre for Australian Weather and Climate Research (CAWCR),
A partnership between the Australian Bureau of Meteorology and CSIRO
email: j.pe...@bom.gov.aumailto:j.pe...@bom.gov.au
phone: +61 (0) 3 9669 4838
fax: +61 (0) 3 9669 4660
mail: GPO Box 1289K, Melbourne, Victoria 3001, Australia
[[alternative HTML version deleted]]
__
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.
[R] 3d topographic map
Hi All- 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? thanks- sherri __ 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.
Re: [R] 3d topographic map
On Jul 25, 2010, at 6:31 PM, sh...@ucar.edu wrote: Hi All- 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? thanks- I suspect most viewers of this message are going to be puzzled. If wireframe and contourplot are not doing it for you, what are the problems? Surely you found references to Lattice and wireframe. The Lattice system has a worked example in Figure 6.11 from this page: http://lmdvr.r-forge.r-project.org/figures/figures.html A search for topographic in one of the r searching facilities is sure to bring multiple hits, . including the usual starting point: https://svn.r-project.org/R/trunk/src/library/graphics/demo/graphics.R So there must be details and issues that you have not chosen to disclose. -- David. __ 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.
