1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
%#
%# fields is a package for analysis of spatial data written for
%# the R software environment.
%# Copyright (C) 2024 Colorado School of Mines
%# 1500 Illinois St., Golden, CO 80401
%# Contact: Douglas Nychka, douglasnychka@gmail.edu,
%#
%# This program is free software; you can redistribute it and/or modify
%# it under the terms of the GNU General Public License as published by
%# the Free Software Foundation; either version 2 of the License, or
%# (at your option) any later version.
%# This program is distributed in the hope that it will be useful,
%# but WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%# GNU General Public License for more details.
%#
%# You should have received a copy of the GNU General Public License
%# along with the R software environment if not, write to the Free Software
%# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
%# or see http://www.r-project.org/Licenses/GPL-2
%##END HEADER
%##END HEADER
\name{glacier}
\alias{glacier}
\docType{data}
\title{Franke's Glacier Elevation Data}
\description{
A moderate size (about 8400 locations) spatial dataset that is well-known in the applied mathematics approximation literature for testing interpolation methods.
}
\usage{data(glacier)}
\format{
The format of \code{glacier} is a list with two components:
\describe{
\item{loc:}{8338x2 matrix of the locations (meters??).}
\item{y:}{A vector of elevations (meters ??).}
}
}
\details{
This data set appears in papers that develop interpolation methods for
scattered data and serves as an interesting bridge to the examples
in applied math that develop radial basis function surface fitting.
The data was originally used by R. Franke.
Unfortunately at this time we can not find any background on where
these data were collected or indeed even the location of this glacier.
However, it is an interesting data set in that it appears that
the elevations are reported along lnes of equal elevation, i.e.
contours, perhaps from a digitization of a
topographic map or survey. It is important to estimate the
surface in a way that the artifacts from discretization are not
present. In the example below the compactly supported kernel
interpolation still has some artifacts.
The glacier data set is available at this website
\url{https://oleg-davydov.de/scat_data.html}
The examples below are useful for comparing different
approximations
to a Gaussian spatial process estimate for the elevation surface.
Of
course in using a stationary covariance ( e.g. the Matern or
Wendland)
these are also radial basis smoothing or interpolation of the
data.
}
\examples{
data( glacier )
# EDA for raw obs:
bubblePlot( glacier$loc, glacier$y, highlight=FALSE, size=.5)
# identifying contour levels. Note this is reported at regular levels
# (Every 25m ???)
table( glacier$y)
# find sigma and rho by maximum likelihood
# for a fixed range
# the default is the Wendland covariance with k=2
# See help(Wendland)
# this takes about 5 minutes
# macbook pro Quad-Core Intel Core i5 8 GB
#options(spam.nearestdistnnz=c(5e7,1e3))
#system.time(
# obj0<- fastTps(glacier$loc, glacier$y,
# theta=2,
# profileLambda=TRUE)
#)
# set.panel(2,2)
# plot( obj0)
# set.panel()
# just evaluate at MLE
# reset default matrix size that the spam pacakge will use.
\dontrun{
options(spam.nearestdistnnz=c(5e7,1e3))
system.time( obj1<-
fastTps(glacier$loc, glacier$y,
theta=2,
lambda= 7.58e-5
)
)
system.time(
look1<- predictSurface( obj1, nx=150, ny=150)
)
imagePlot( look1)
system.time(
out<- simLocal.spatialProcess(obj1, M=3, nx=150, ny=150)
)
set.panel( 2,2)
imagePlot( look1)
zlim<- range( out$z, na.rm=TRUE)
for( k in 1:3){
imagePlot(out$x, out$y, out$z[,,k], zlim=zlim)
}
# near interpolation surface using Matern smoothness .5
system.time(
obj2<- spatialProcess(glacier$loc, glacier$y,
aRange = 1.5,
lambda = 1e-5,
smoothness = .5)
)
system.time(
out<- simLocal.spatialProcess(obj2, M=3, nx=150, ny=150,
fast=TRUE)
)
set.panel( 2,2)
imagePlot( look1)
zlim<- range( out$z, na.rm=TRUE)
for( k in 1:3){
imagePlot(out$x, out$y, out$z[,,k], zlim=zlim)
}
% test out fast predict algorithm verses exact
% note speedup of about 15 times
system.time(
look2<- predictSurface.mKrig( obj2, nx=150, ny=150,
fast=TRUE, NNSize=5)
)
system.time(
look2B<- predictSurface( obj2, nx=150, ny=150,
fast=FALSE)
)
err<- c((look2$z - look2B$z)/look2B$z)
stats( log10( abs(err) ) )
# some error plots ( percent relative error)
imagePlot(look2$x, look2$y, 100*(look2$z - look2B$z)/look2B$z )
imagePlot(look2$x, look2$y, 100*(look1$z - look2B$z)/look2B$z )
} % end do not run
} % end examples
\keyword{datasets}
|
