|Contact Us | Visit Us | UCAR People Search||Numerics| Assimilation| Turbulence| Statistics|
For all of these examples you should load the fields package and also the datasets. (You only need to do this once when starting R.) Fields is available as a contributed package in the Comprehensive R Archive Network CRAN
Note that the perspective and image plots actually use the lon lat coordinates. But the Euclidean transformation is used for surface fitting and estimating derivatives.
# find spacing in meters of the lon and lat grid.
# and scale of longitude coordinates to give equal distance
# Because this region is so small we will ignore the curvature of
# coordinates on sphere.
dlat<- RifleDEM$y- RifleDEM$y
dx<- rdist.earth( cbind( RifleDEM$x[1:2], RifleDEM$y[c(1,1)] ),
dy<- rdist.earth( cbind( RifleDEM$x[c(1,1)], RifleDEM$y[1:2] ) ,
dlon<- RifleDEM$x- RifleDEM$x
dlat<- RifleDEM$y- RifleDEM$y
# Create all the grid locations explicitly
# converting to a approximate cartesian grid in meters
make.surface.grid( list( x=(RifleDEM$x- RifleDEM$x)*dx/dlon,
y= (RifleDEM$y- RifleDEM$y)*dy/dlat) )-> Rgrid
# Rtrial is the ski run in the meter Euclidean coordinates.
(rifle.trail[,1] - RifleDEM$x)*dx/dlon,
(rifle.trail[,2] - RifleDEM$y)*dy/dlat )
# to find slopes need to parameterthe 1-d curve of the ski run
# by arc length
N<- nrow( rifle.trail)
rD<- sqrt(diff( Rtrail[,1])**2 + diff( Rtrail[,2])**2)
trail.D<- c(0,cumsum( rD))
#First remove spatial linear drift by least squares fit before finding
# the variogram.
zR<- lsfit( Rgrid, c( RifleDEM$z))$residuals
zR<- matrix( zR, ncol= ncol( RifleDEM$z), nrow= nrow( RifleDEM$z))
# see help(vgram.matrix) for details on this function
# R is set to look at all pairs within a distance of about 6.5 grid spacings
vgram.matrix( zR, R= 65, dx=dx, dy=dy) -> look
# A nice plot converting distance in scaled lon/lat to meters.
plot( look$d.full, look$vgram.full, xlab="Separation Distance (m)",
ylab="Variogram", xlim=c(0, 35),ylim=c(0,32), col=2)
# Kriging with compactly supported Wendland covariance function
# The range (theta) controls support.
mKrig( Rgrid, c( RifleDEM$z), cov.function="wendland.cov", k=2,
lambda=lambda,mean.neighbor=50 )-> rifle.fit
rifle.trail.elev<- predict(rifle.fit, Rtrail)
# computing slope based on 1-d space curveIt can be checked that the two estimates of slope are very similar.
# and interpolating cubic spline
trail.slope<- splint( trail.D,rifle.trail.elev, trail.D, derivative=1)
trail.degrees<- - atan( trail.slope)*(360/(2*pi))
# The correct way to find by slope taking partials of elevation surface
rifle.trail.der<- predict(rifle.fit, Rtrail,derivative=1)
# First find unit tangent vector along run
xp<- splint( trail.D,rifle.trail[,1], trail.D, derivative=1)
yp<- splint( trail.D,rifle.trail[,2], trail.D, derivative=1)
R<- sqrt( xp**2 + yp**2)
# find directional derivative based on tangent vector
trail.slope2<- xp*rifle.trail.der[,1] + yp*rifle.trail.der[,2]
trail.degrees2<- - atan( trail.slope2)*(360/(2*pi))
matplot( trail.D, cbind(trail.degrees,trail.degrees2),
type="l", xlab="Distance (m)",
ylab="Slope degrees",lwd=2, col=c(2,3),lty=1 )
# Just for fun ... add in elevation drop
par( usr=c(0,2000,2800,3500) )
lines( trail.D, rifle.trail.elev,type="l", xlab="Distance (m)",
ylab="Elevation (m)",col=4, lwd=1.5, lty=2)