Revision ae79872e
Added by Jim Regetz over 12 years ago
- ID ae79872e2385a09773e7f2f3c8e45133190beffd
| terrain/research/oregon/arcgis/v2/multiscalesmooth9a_clean.R | ||
|---|---|---|
| 67 | 67 |
# set aggregation factor |
| 68 | 68 |
AGG.FACTOR <- 3 |
| 69 | 69 |
|
| 70 |
# naming: |
|
| 71 |
# h - the value being smoothed/interpolated |
|
| 72 |
# vg - total variance of group of data, or of individual measurement |
|
| 73 |
# v_bg - variance between groups |
|
| 74 |
# v_wg - variance within groups |
|
| 75 |
# wa - weighting for aggregation, based on total variance |
|
| 76 |
# vm - variance of the calculated mean |
|
| 77 |
# mv - mean of finer scale variances |
|
| 78 |
# n - effective number of measurements |
|
| 79 |
|
|
| 80 |
# NB - only calculating sample variances here, not variances of estimated means. |
|
| 81 |
# Also note that v0_bg is an uncertainty, not a sample variance |
|
| 82 |
# and v1_bg is total variances, but both are labelled as "between-group" to simplify the smoothing |
|
| 83 |
|
|
| 84 | 70 |
# expand grid to accommodate integer number of cells at coarsest |
| 85 | 71 |
# level, by adding roungly equal number of cells on either side |
| 86 |
# (with on extra on top/right if an odd number of cells needs to be |
|
| 72 |
# (with one extra on top/right if an odd number of cells needs to be
|
|
| 87 | 73 |
# added) |
| 88 | 74 |
max.size <- AGG.FACTOR^NUM.LEVELS |
| 89 | 75 |
addx <- if (0 < (extra<-ncol(ingrid)%%max.size)) {
|
| ... | ... | |
| 96 | 82 |
} else {
|
| 97 | 83 |
0 |
| 98 | 84 |
} |
| 99 |
h <- list(expand(ingrid, extent(
|
|
| 85 |
z <- list(expand(ingrid, extent(
|
|
| 100 | 86 |
xmin(ingrid)-floor(addx)*xres(ingrid), |
| 101 | 87 |
xmax(ingrid)+ceiling(addx)*xres(ingrid), |
| 102 | 88 |
ymin(ingrid)-floor(addy)*yres(ingrid), |
| 103 | 89 |
ymax(ingrid)+ceiling(addy)*yres(ingrid) |
| 104 | 90 |
))) |
| 105 |
v <- list(calc(h[[1]], function(x) ifelse(!is.na(x), sd^2, NA))) |
|
| 91 |
v <- list(calc(z[[1]], function(x) ifelse(!is.na(x), sd^2, NA))) |
|
| 92 |
|
|
| 93 |
# create grids |
|
| 94 |
|
|
| 95 |
# NB - only calculating sample variances here, not variances of estimated means. |
|
| 96 |
# Also note that v0_bg is an uncertainty, not a sample variance |
|
| 97 |
# and v1_bg is total variances, but both are labelled as "between-group" to simplify the smoothing |
|
| 106 | 98 |
|
| 107 |
vg = v[[1]] |
|
| 99 |
# set initial "group variance" to individual msmt variance (noise) |
|
| 100 |
v.g = v[[1]] |
|
| 101 |
# weighting for aggregation, based on total variance |
|
| 108 | 102 |
w <- calc(v[[1]], function(x) ifelse(is.na(x), 0, 1/x)) |
| 109 | 103 |
wsq = w^2 |
| 110 |
n <- calc(h[[1]], function(x) ifelse(!is.na(x), 1, 0)) |
|
| 104 |
# effective number of measurements |
|
| 105 |
n <- calc(z[[1]], function(x) ifelse(!is.na(x), 1, 0)) |
|
| 111 | 106 |
|
| 112 | 107 |
bigvar <- cellStats(v[[1]], max) |
| 113 | 108 |
|
| ... | ... | |
| 118 | 113 |
|
| 119 | 114 |
message("Aggregate from ", i-1, " to ", i)
|
| 120 | 115 |
|
| 121 |
vg.prev <- vg |
|
| 116 |
# make copies of previous (finer scale) grids |
|
| 117 |
v.g.prev <- v.g |
|
| 122 | 118 |
w.prev <- w |
| 123 | 119 |
wsq.prev <- wsq |
| 124 | 120 |
n.prev <- n |
| 125 | 121 |
|
| 122 |
# calc neighborhood weights, num cells, effective num cells |
|
| 126 | 123 |
w <- aggregate(w.prev, 3, sum) |
| 127 | 124 |
wsq <- aggregate(wsq.prev, 3, sum) |
| 128 | 125 |
n <- aggregate(n.prev, 3, sum) |
| 129 |
neff <- w^2 / wsq |
|
| 130 |
h[[i]] <- aggregate(w.prev * h[[i-1]], 3, sum) / w |
|
| 131 |
hdiff <- h[[i-1]] - disaggregate(h[[i]], 3) |
|
| 132 |
vbg <- aggregate(w.prev * hdiff^2, 3, sum) / w |
|
| 126 |
n.eff <- w^2 / wsq |
|
| 127 |
# calc variance-weighted neighborhood mean |
|
| 128 |
z[[i]] <- aggregate(w.prev * z[[i-1]], 3, sum) / w |
|
| 129 |
# calc between-cell variance, taken over neighborhood |
|
| 130 |
hdiff <- z[[i-1]] - disaggregate(z[[i]], 3) |
|
| 131 |
v.bg <- aggregate(w.prev * hdiff^2, 3, sum) / w |
|
| 132 |
# calc wtd avg of within-cell variance, taken over neighborhood |
|
| 133 | 133 |
if (i==2) {
|
| 134 |
vwg <- n - n # zero, but with window and cell size set for us |
|
| 134 |
v.wg <- n - n # zero, but with window and cell size set for us
|
|
| 135 | 135 |
} else {
|
| 136 |
vwg <- aggregate(w.prev * vg.prev, 3, sum) / w
|
|
| 136 |
v.wg <- aggregate(w.prev * v.g.prev, 3, sum) / w
|
|
| 137 | 137 |
} |
| 138 |
vg <- vbg + vwg |
|
| 139 |
vm <- 1 / w |
|
| 138 |
# calc total group variance |
|
| 139 |
# ~ total variance of cell values in the underlying neighborhood |
|
| 140 |
v.g <- v.bg + v.wg |
|
| 141 |
# calc variance of the mean for the neighborhood |
|
| 142 |
v.m <- 1 / w |
|
| 143 |
# calc mean noise variance (mean of finer scale variances) |
|
| 140 | 144 |
mv <- n / w |
| 141 | 145 |
|
| 142 |
chisq <- 1 + chisqa / sqrt(neff - 1) + chisqb / (neff - 1) |
|
| 143 |
v[[i]] <- overlay(vg, chisq, mv, vm, |
|
| 144 |
fun=function(vg, chisq, mv, vm) ifelse(vg/chisq < mv, vm, vg)) |
|
| 146 |
chisq <- 1 + chisqa / sqrt(n.eff - 1) + chisqb / (n.eff - 1) |
|
| 147 |
# set coarsened cell variances: if group variance is small |
|
| 148 |
# relative to noise variance, use variance of the mean instead |
|
| 149 |
# of group variance |
|
| 150 |
v[[i]] <- overlay(v.g, chisq, mv, v.m, |
|
| 151 |
fun=function(v.g, chisq, mv, v.m) ifelse(v.g/chisq < mv, v.m, v.g)) |
|
| 145 | 152 |
|
| 146 | 153 |
} |
| 147 | 154 |
|
| 148 | 155 |
bigvar <- bigvar * 10 |
| 149 | 156 |
|
| 150 |
#copy h[NUM.LEVELS] hs[NUM.LEVELS]
|
|
| 157 |
#copy z[NUM.LEVELS] hs[NUM.LEVELS]
|
|
| 151 | 158 |
#copy v[NUM.LEVELS] vs[NUM.LEVELS] |
| 152 |
#kill h[NUM.LEVELS]
|
|
| 159 |
#kill z[NUM.LEVELS]
|
|
| 153 | 160 |
#kill v[NUM.LEVELS] |
| 154 | 161 |
#setcell hs[NUM.LEVELS] |
| 155 | 162 |
#setwindow hs[NUM.LEVELS] |
| 156 | 163 |
|
| 157 | 164 |
# smooth, refine and combine each layer in turn |
| 158 |
#hs <- h[[NUM.LEVELS+1]]
|
|
| 165 |
#hs <- z[[NUM.LEVELS+1]]
|
|
| 159 | 166 |
#vs <- v[[NUM.LEVELS+1]] |
| 160 |
hs <- h
|
|
| 167 |
hs <- z
|
|
| 161 | 168 |
vs <- v |
| 162 | 169 |
#todo: is this the same as circle with radius 2 in arc??? |
| 163 | 170 |
#circle2 <- matrix(c(0,1,0, 1,1,1, 0,1,0), nrow=3) |
| ... | ... | |
| 169 | 176 |
message("Refine from ", j, " to ", j-1)
|
| 170 | 177 |
|
| 171 | 178 |
# for the first stage where the coarser grid is refined and smoothed, set window to the coarse grid |
| 172 |
#setcell h[j-1]
|
|
| 179 |
#setcell z[j-1]
|
|
| 173 | 180 |
#setwindow maxof |
| 174 | 181 |
|
| 175 |
# create smoothed higher resolution versions of h and v_bg, hopefully with no nulls!
|
|
| 182 |
# create smoothed higher resolution versions of z and v_bg, hopefully with no nulls!
|
|
| 176 | 183 |
# [suppressing warnings to avoid .couldBeLonLat complaints] |
| 177 | 184 |
suppressWarnings({
|
| 178 | 185 |
hs.tmp <- disaggregate(focal(hs[[j]], w=circle2, fun=mean, |
| ... | ... | |
| 181 | 188 |
pad=TRUE, padValue=NA, na.rm=TRUE), 3) |
| 182 | 189 |
}) |
| 183 | 190 |
|
| 184 |
# create no-null version of finer h and v
|
|
| 185 |
h_c <- calc(h[[j-1]], function(x) ifelse(is.na(x), 0, x))
|
|
| 191 |
# create no-null version of finer z and v
|
|
| 192 |
h_c <- calc(z[[j-1]], function(x) ifelse(is.na(x), 0, x))
|
|
| 186 | 193 |
v_c <- calc(v[[j-1]], function(x) ifelse(is.na(x), bigvar, x)) |
| 187 | 194 |
|
| 188 | 195 |
# combine two values using least variance |
Also available in: Unified diff
tweaked var names to match Gallant pub; added comments