1
|
# Multiscale raster smoother function designed for use with noisy DEMs.
|
2
|
# Original algorithm and AML implementation by John Gallant. Translation
|
3
|
# to R (relying on raster package functionality) by Jim Regetz, based on
|
4
|
# AML script 'multiscalesmooth9a_clean.aml'.
|
5
|
#
|
6
|
# Jim Regetz
|
7
|
# NCEAS
|
8
|
# Created on 22-Mar-2012
|
9
|
|
10
|
library(raster)
|
11
|
|
12
|
multiscalesmooth <- function(ingrid, sd=0.0001, alpha=0.05, bbox) {
|
13
|
# ingrid: RasterLayer to smooth
|
14
|
# sd: numeric constant or RasterLayer
|
15
|
# alpha: alpha level for chi-square critical value
|
16
|
# bbox: optional extent object used to subset ingrid
|
17
|
|
18
|
# subset ingrid if desired
|
19
|
if (!missing(bbox)) {
|
20
|
ingrid <- crop(ingrid, bbox)
|
21
|
}
|
22
|
|
23
|
# set number of aggregation levels and neighborhood size
|
24
|
NUM.LEVELS <- 4
|
25
|
NUM.CELLS <- 3
|
26
|
|
27
|
# expand grid to accommodate integer number of cells at coarsest
|
28
|
# level, by adding roungly equal number of cells on either side
|
29
|
# (with one extra on top/right if an odd number of cells is needed)
|
30
|
max.size <- NUM.CELLS^NUM.LEVELS
|
31
|
addx <- if (0 < (extra<-ncol(ingrid)%%max.size)) {
|
32
|
(max.size-extra)/2
|
33
|
} else {
|
34
|
0
|
35
|
}
|
36
|
addy <- if (0 < (extra<-nrow(ingrid)%%max.size)) {
|
37
|
(max.size-extra)/2
|
38
|
} else {
|
39
|
0
|
40
|
}
|
41
|
full.extent <- extent(
|
42
|
xmin(ingrid) - floor(addx) * xres(ingrid),
|
43
|
xmax(ingrid) + ceiling(addx) * xres(ingrid),
|
44
|
ymin(ingrid) - floor(addy) * yres(ingrid),
|
45
|
ymax(ingrid) + ceiling(addy) * yres(ingrid)
|
46
|
)
|
47
|
|
48
|
# create grids
|
49
|
|
50
|
# create lists to hold the series of successively coarsened grids of
|
51
|
# values and variances, respectively; the first element of each is
|
52
|
# the grid at the original input resolution
|
53
|
# ...insert initial grid of values, but expanded to full extent
|
54
|
z <- list(expand(ingrid, full.extent))
|
55
|
# ...insert initial grid of variances
|
56
|
if (is.numeric(sd) && length(sd)==1) {
|
57
|
v <- list(calc(z[[1]], function(x) ifelse(!is.na(x), sd^2, NA)))
|
58
|
} else if (class(sd)=="RasterLayer") {
|
59
|
if (identical(extent(sd), extent(ingrid))) {
|
60
|
v <- list(overlay(z[[1]], sd, fun=function(z, sd)
|
61
|
ifelse(!is.na(z), sd^2, NA)))
|
62
|
} else {
|
63
|
stop("sd raster extent differs from ingrid extent")
|
64
|
}
|
65
|
} else {
|
66
|
stop("sd must be a single number or a RasterLayer")
|
67
|
}
|
68
|
|
69
|
# set initial "group variance" to individual msmt variance (noise)
|
70
|
v.g = v[[1]]
|
71
|
# weights for aggregation, based on total variance
|
72
|
w <- calc(v[[1]], function(x) ifelse(is.na(x), 0, 1/x))
|
73
|
# squared weights
|
74
|
wsq = w^2
|
75
|
# effective number of measurements
|
76
|
n <- calc(z[[1]], function(x) ifelse(!is.na(x), 1, 0))
|
77
|
|
78
|
# aggregate to broader scales
|
79
|
for (i in 1+seq.int(NUM.LEVELS)) {
|
80
|
|
81
|
message("Aggregate from ", i-1, " to ", i)
|
82
|
|
83
|
# make copies of previous (finer scale) grids
|
84
|
v.g.prev <- v.g
|
85
|
w.prev <- w
|
86
|
wsq.prev <- wsq
|
87
|
n.prev <- n
|
88
|
|
89
|
# calc neighborhood weights, num cells, effective num cells
|
90
|
w <- aggregate(w.prev, 3, sum)
|
91
|
wsq <- aggregate(wsq.prev, 3, sum)
|
92
|
n <- aggregate(n.prev, 3, sum)
|
93
|
n.eff <- w^2 / wsq
|
94
|
# calc variance-weighted neighborhood mean
|
95
|
z[[i]] <- aggregate(w.prev * z[[i-1]], 3, sum) / w
|
96
|
# calc between-cell variance, taken over neighborhood
|
97
|
zdiff <- z[[i-1]] - disaggregate(z[[i]], 3)
|
98
|
v.bg <- aggregate(w.prev * zdiff^2, 3, sum) / w
|
99
|
# calc wtd avg of within-cell variance, taken over neighborhood
|
100
|
if (i==2) {
|
101
|
v.wg <- n - n # zero, but with correct window and cell size
|
102
|
} else {
|
103
|
v.wg <- aggregate(w.prev * v.g.prev, 3, sum) / w
|
104
|
}
|
105
|
# calc total group variance
|
106
|
# ~ total variance of cell values in the underlying neighborhood
|
107
|
v.g <- v.bg + v.wg
|
108
|
# calc variance of the mean for the neighborhood
|
109
|
v.m <- 1 / w
|
110
|
# calc mean noise variance (mean of finer scale variances)
|
111
|
mv <- n / w
|
112
|
|
113
|
# calc chisq critical values
|
114
|
chisq <- calc(n.eff, function(n) qchisq(alpha, n-1,
|
115
|
lower=FALSE)/(n-1))
|
116
|
# set coarsened cell variances: if group variance is small
|
117
|
# relative to noise variance, use variance of the mean instead
|
118
|
# of group variance
|
119
|
v[[i]] <- overlay(v.g, chisq, mv, v.m,
|
120
|
fun=function(v.g, chisq, mv, v.m) ifelse(v.g/chisq < mv, v.m, v.g))
|
121
|
|
122
|
}
|
123
|
|
124
|
# get arbitrarily large value to fill null variance cells
|
125
|
bigvar <- cellStats(v[[1]], max) * 10
|
126
|
|
127
|
# prep for refinement phase
|
128
|
z.smooth <- z[[NUM.LEVELS+1]]
|
129
|
v.smooth <- v[[NUM.LEVELS+1]]
|
130
|
# create weight matrix equivalent to ArcGIS circle with radius 2
|
131
|
circle2 <- matrix(c(0,0,1,0,0,
|
132
|
0,1,1,1,0,
|
133
|
1,1,1,1,1,
|
134
|
0,1,1,1,0,
|
135
|
0,0,1,0,0), nrow=5)
|
136
|
circle2[circle2==0] <- NA
|
137
|
|
138
|
# refine, smooth, and combine each layer in turn
|
139
|
for (j in 1+rev(seq.int(NUM.LEVELS))) {
|
140
|
|
141
|
message("Refining from ", j, " to ", j-1)
|
142
|
|
143
|
# create smoothed higher resolution versions of z and v
|
144
|
zs <- focal(disaggregate(z.smooth, 3), w=circle2, fun=mean,
|
145
|
pad=TRUE, padValue=NA, na.rm=TRUE)
|
146
|
vs <- focal(disaggregate(v.smooth, 3), w=circle2, fun=mean,
|
147
|
pad=TRUE, padValue=NA, na.rm=TRUE)
|
148
|
|
149
|
# create no-null version of finer z and v
|
150
|
z_c <- calc(z[[j-1]], function(x) ifelse(is.na(x), 0, x))
|
151
|
v_c <- calc(v[[j-1]], function(x) ifelse(is.na(x), bigvar, x))
|
152
|
|
153
|
# explicitly clean up, in case it helps
|
154
|
z[[j-1]] <- NULL
|
155
|
v[[j-1]] <- NULL
|
156
|
|
157
|
# combine two values using least variance
|
158
|
z.smooth <- (z_c/v_c + zs/vs ) / (1/v_c + 1/vs)
|
159
|
v.smooth <- 1 / (1/v_c + 1/vs)
|
160
|
|
161
|
}
|
162
|
|
163
|
result <- crop(stack(z.smooth, v.smooth, extent(ingrid))
|
164
|
layerNames(result) <- c("zs", "vs")
|
165
|
return(result)
|
166
|
|
167
|
}
|