Splines on the sphere
gam can use isotropic smooths on the sphere, via terms like
s(la,lo,bs="sos",m=2,k=100). There must be exactly 2 arguments to such a smooth.
The first is taken to be latitude (in degrees) and the second longitude (in degrees).
m (default 0) is an integer in the range -1 to 4 determining the order of the penalty used.
For m>0, (m+2)/2 is the penalty order, with m=2 equivalent to the usual second
derivative penalty. m=0 signals to use the 2nd order spline on the sphere, computed by
Wendelberger's (1981) method. m = -1 results in a Duchon.spline being
used (with m=2 and s=1/2), following an unpublished suggestion of Jean Duchon.
k (default 50) is the basis dimension.
## S3 method for class 'sos.smooth.spec' smooth.construct(object, data, knots) ## S3 method for class 'sos.smooth' Predict.matrix(object, data)
object |
a smooth specification object, usually generated by a term |
data |
a list containing just the data (including any |
knots |
a list containing any knots supplied for basis setup — in same order and with same names as |
For m>0, the smooths implemented here are based on the pseudosplines on the sphere of Wahba (1981)
(there is a correction of table 1 in 1982, but the correction has a misprint in the definition of A — the A
given in the 1981 paper is correct). For m=0 (default) then a second order spline on the sphere is used which is the
analogue of a second order thin plate spline in 2D: the computation is based on Chapter 4 of Wendelberger, 1981.
Optimal low rank approximations are obtained using exactly the approach given in Wood (2003). For m = -1 a smooth of the general type discussed in Duchon (1977) is used: the sphere is embedded in a 3D Euclidean space, but
smoothing employs a penalty based on second derivatives (so that locally as the smoothing parameter tends
to zero we recover a "normal" thin plate spline on the tangent space). This is an unpublished suggestion of Jean
Duchon. m = -2 is the same but with first derivative penalization.
Note that the null space of the penalty is always the space of constant functions on the sphere, whatever the order of penalty.
This class has a plot method, with 3 schemes. scheme==0 plots one hemisphere of the sphere, projected onto a circle.
The plotting sphere has the north pole at the top, and the 0 meridian running down the middle of the plot, and towards
the viewer. The smoothing sphere is rotated within the plotting sphere, by specifying the location of its pole in the
co-ordinates of the viewing sphere. theta, phi give the longitude and latitude of the smoothing sphere pole
within the plotting sphere (in plotting sphere co-ordinates). (You can visualize the smoothing sphere as a globe, free to
rotate within the fixed transparent plotting sphere.) The value of the smooth is shown by a heat map overlaid with a
contour plot. lat, lon gridlines are also plotted.
scheme==1 is as scheme==0, but in black and white, without the image plot. scheme>1 calls the default
plotting method with scheme decremented by 2.
An object of class "sos.smooth". In addition to the usual elements of a
smooth class documented under smooth.construct, this object will contain:
Xu |
A matrix of the unique covariate combinations for this smooth (the basis is constructed by first stripping out duplicate locations). |
UZ |
The matrix mapping the parameters of the reduced rank spline back to the parameters of a full spline. |
Simon Wood simon.wood@r-project.org, with help from Grace Wahba (m=0 case) and Jean Duchon (m = -1 case).
Wahba, G. (1981) Spline interpolation and smoothing on the sphere. SIAM J. Sci. Stat. Comput. 2(1):5-16
Wahba, G. (1982) Erratum. SIAM J. Sci. Stat. Comput. 3(3):385-386.
Wendelberger, J. (1981) PhD Thesis, University of Winsconsin.
Wood, S.N. (2003) Thin plate regression splines. J.R.Statist.Soc.B 65(1):95-114
require(mgcv)
set.seed(0)
n <- 400
f <- function(la,lo) { ## a test function...
sin(lo)*cos(la-.3)
}
## generate with uniform density on sphere...
lo <- runif(n)*2*pi-pi ## longitude
la <- runif(3*n)*pi-pi/2
ind <- runif(3*n)<=cos(la)
la <- la[ind];
la <- la[1:n]
ff <- f(la,lo)
y <- ff + rnorm(n)*.2 ## test data
## generate data for plotting truth...
lam <- seq(-pi/2,pi/2,length=30)
lom <- seq(-pi,pi,length=60)
gr <- expand.grid(la=lam,lo=lom)
fz <- f(gr$la,gr$lo)
zm <- matrix(fz,30,60)
require(mgcv)
dat <- data.frame(la = la *180/pi,lo = lo *180/pi,y=y)
## fit spline on sphere model...
bp <- gam(y~s(la,lo,bs="sos",k=60),data=dat)
## pure knot based alternative...
ind <- sample(1:n,100)
bk <- gam(y~s(la,lo,bs="sos",k=60),
knots=list(la=dat$la[ind],lo=dat$lo[ind]),data=dat)
b <- bp
cor(fitted(b),ff)
## plot results and truth...
pd <- data.frame(la=gr$la*180/pi,lo=gr$lo*180/pi)
fv <- matrix(predict(b,pd),30,60)
par(mfrow=c(2,2),mar=c(4,4,1,1))
contour(lom,lam,t(zm))
contour(lom,lam,t(fv))
plot(bp,rug=FALSE)
plot(bp,scheme=1,theta=-30,phi=20,pch=19,cex=.5)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.