1 |
86eb951b
|
Jim Regetz
|
/* run in GRID, not ARC
|
2 |
|
|
/* &run .aml ingrid sd prob bbox
|
3 |
|
|
/* sd 0.0001 prob 0.05
|
4 |
|
|
|
5 |
|
|
&args ingrid sd prob bbox:rest
|
6 |
|
|
|
7 |
|
|
/* 9a - limit to 4 steps, see if that has any significant deterioration of smoothing performance. Should fix
|
8 |
|
|
/* the problem with islands and headlands - turns out also need to remove water except
|
9 |
|
|
/* for one cell adjacent to land, and give that a higher uncertainty. See cluster_multiscalesmooth9a_clean_216
|
10 |
|
|
|
11 |
|
|
/* Version 9:
|
12 |
|
|
/* Focus on implementing identical algorithm to directsmooth2 using multiscale method i.e. aggregating by factor of 3
|
13 |
|
|
/* from already aggregated data, rather than from original resolution each time.
|
14 |
|
|
|
15 |
|
|
/* Version 8:
|
16 |
|
|
/* WARNING: have not yet checked that the additional weighting of the gaussian smoothing is not messing with
|
17 |
|
|
/* the calculations of variance etc.
|
18 |
|
|
/* Replaced simple 3x3 aggregation with gaussian smoothing
|
19 |
|
|
/* Kernel is chosen to give appropriate level of smoothing and produce a fairly accurate approximation of the smoothed
|
20 |
|
|
/* surface by interpolation of the coarser scale smoothed values
|
21 |
|
|
/* Details in gaussian.xls on john's laptop
|
22 |
|
|
|
23 |
|
|
/* Version 7:
|
24 |
|
|
/* Further reworking of how the final values are selected - a mean is a candidate if its associated variance
|
25 |
|
|
/* is less than the mean sample uncertainty, and the mean with the lowest variance among the candidates is the chosen one.
|
26 |
|
|
/* Implement that in the nested sense by taking the lowest group variance divided by the chi^2 value, and its associated mean variance,
|
27 |
|
|
/* and if that is lower than the data point variance the
|
28 |
|
|
|
29 |
|
|
/* approximate critical value of chi^2/N with N degrees of freedom at 5% level as 1 + 2.45/sqrt(N) + 0.55/N
|
30 |
|
|
/* for 1% level use 1 + 3.4/sqrt(N) + 2.4/N
|
31 |
|
|
|
32 |
|
|
/* Version 6:
|
33 |
|
|
/* Done from scratch after careful working through of theory.
|
34 |
|
|
/* ingrid is the (potentially sparse) grid of data to be smoothed and
|
35 |
|
|
/* interpolated, which can be of different extent to the output
|
36 |
|
|
/* (resolution is assumed to be the same, could adapt this to relax that)
|
37 |
|
|
/*
|
38 |
|
|
/* var can be a constant or a grid
|
39 |
|
|
/*
|
40 |
|
|
/* bbox can be either a grid name or the 'xmin ymin xmax ymax' parameters
|
41 |
|
|
/* for setwindow
|
42 |
|
|
|
43 |
|
|
&type NB - using standard deviation as noise specification now, not variance!
|
44 |
|
|
|
45 |
|
|
|
46 |
|
|
/* set up chisq parameters
|
47 |
|
|
&sv chisqa = [calc 2.807 - 0.6422 * [log10 %prob% ] - 3.410 * %prob% ** 0.3411 ]
|
48 |
|
|
&sv chisqb = [calc -5.871 - 3.675 * [log10 %prob% ] + 4.690 * %prob% ** 0.3377 ]
|
49 |
|
|
&type chisq parameters %chisqa% %chisqb%
|
50 |
|
|
|
51 |
|
|
|
52 |
|
|
setcell %ingrid%
|
53 |
|
|
|
54 |
|
|
/* work out maximum of ingrid and bbox extents
|
55 |
|
|
setwindow [unquote %bbox%] %ingrid%
|
56 |
|
|
bboxgrid = 1
|
57 |
|
|
setwindow maxof
|
58 |
|
|
workgrid = %ingrid% + bboxgrid
|
59 |
|
|
setwindow workgrid
|
60 |
|
|
kill workgrid
|
61 |
|
|
|
62 |
|
|
/* naming:
|
63 |
|
|
/* h - the value being smoothed/interpolated
|
64 |
|
|
/* vg - total variance of group of data, or of individual measurement
|
65 |
|
|
/* v_bg - variance between groups
|
66 |
|
|
/* v_wg - variance within groups
|
67 |
|
|
/* wa - weighting for aggregation, based on total variance
|
68 |
|
|
/* vm - variance of the calculated mean
|
69 |
|
|
/* mv - mean of finer scale variances
|
70 |
|
|
/* n - effective number of measurements
|
71 |
|
|
|
72 |
|
|
|
73 |
|
|
/* NB - only calculating sample variances here, not variances of estimated means.
|
74 |
|
|
/* Also note that v0_bg is an uncertainty, not a sample variance
|
75 |
|
|
/* and v1_bg is total variances, but both are labelled as "between-group" to simplify the smoothing
|
76 |
|
|
|
77 |
|
|
h0 = %ingrid%
|
78 |
|
|
v0 = con(^ isnull(h0), sqr(%sd%))
|
79 |
|
|
vg0 = v0
|
80 |
|
|
w0 = con(isnull(v0), 0, 1.0 / v0)
|
81 |
|
|
wsq0 = sqr(w0)
|
82 |
|
|
n0 = con(^ isnull(h0), 1, 0)
|
83 |
|
|
|
84 |
|
|
&describe v0
|
85 |
|
|
&sv bigvar %grd$zmax%
|
86 |
|
|
|
87 |
|
|
setcell minof
|
88 |
|
|
|
89 |
|
|
/* aggregate to broader scales
|
90 |
|
|
&sv i 1
|
91 |
|
|
&sv done .false.
|
92 |
|
|
&describe h0
|
93 |
|
|
|
94 |
|
|
&do &until %done%
|
95 |
|
|
&sv j [calc %i% - 1]
|
96 |
|
|
|
97 |
|
|
&type Aggregate from %j% to %i%
|
98 |
|
|
|
99 |
|
|
&describe h%j%
|
100 |
|
|
&sv cell3 [calc %grd$dx% * 3]
|
101 |
|
|
&describe h0
|
102 |
|
|
&sv nx0 [round [calc %grd$xmin% / %cell3% - 0.5]]
|
103 |
|
|
&sv ny0 [round [calc %grd$ymin% / %cell3% - 0.5]]
|
104 |
|
|
&sv nx1 [round [calc %grd$xmax% / %cell3% + 0.5]]
|
105 |
|
|
&sv ny1 [round [calc %grd$ymax% / %cell3% + 0.5]]
|
106 |
|
|
&sv x0 [calc ( %nx0% - 0.5 ) * %cell3%]
|
107 |
|
|
&sv y0 [calc ( %ny0% - 0.5 ) * %cell3%]
|
108 |
|
|
&sv x1 [calc ( %nx1% + 0.5 ) * %cell3%]
|
109 |
|
|
&sv y1 [calc ( %ny1% + 0.5 ) * %cell3%]
|
110 |
|
|
setwindow %x0% %y0% %x1% %y1%
|
111 |
|
|
|
112 |
|
|
w%i% = aggregate(w%j%, 3, sum)
|
113 |
|
|
wsq%i% = aggregate(wsq%j%, 3, sum)
|
114 |
|
|
n%i% = aggregate(n%j%, 3, sum)
|
115 |
|
|
neff%i% = w%i% * w%i% / wsq%i%
|
116 |
|
|
h%i% = aggregate(w%j% * h%j%, 3, sum) / w%i%
|
117 |
|
|
vbg%i% = aggregate(w%j% * sqr(h%j% - h%i%), 3, sum) / w%i%
|
118 |
|
|
&if %i% eq 1 &then vwg%i% = n%i% - n%i% /* zero, but with window and cell size set for us
|
119 |
|
|
&else vwg%i% = aggregate(w%j% * vg%j%, 3, sum) / w%i%
|
120 |
|
|
vg%i% = vbg%i% + vwg%i%
|
121 |
|
|
vm%i% = 1.0 / w%i%
|
122 |
|
|
mv%i% = n%i% / w%i%
|
123 |
|
|
|
124 |
|
|
chisq%i% = 1 + %chisqa% / sqrt(neff%i% - 1) + %chisqb% / (neff%i% - 1)
|
125 |
|
|
v%i% = con(vg%i% / chisq%i% < mv%i%, vm%i%, vg%i%)
|
126 |
|
|
|
127 |
|
|
/* remove everything except h%i% and v%i%
|
128 |
|
|
kill w%j%
|
129 |
|
|
kill wsq%j%
|
130 |
|
|
kill n%j%
|
131 |
|
|
kill neff%i%
|
132 |
|
|
kill vbg%i%
|
133 |
|
|
kill vwg%i%
|
134 |
|
|
kill vg%j%
|
135 |
|
|
kill vm%i%
|
136 |
|
|
kill mv%i%
|
137 |
|
|
kill chisq%i%
|
138 |
|
|
|
139 |
|
|
&sv done %i% eq 4
|
140 |
|
|
|
141 |
|
|
&sv i [calc %i% + 1]
|
142 |
|
|
&end
|
143 |
|
|
|
144 |
|
|
|
145 |
|
|
&sv maxstep [calc %i% - 1]
|
146 |
|
|
&sv bigvar [calc %bigvar% * 10]
|
147 |
|
|
|
148 |
|
|
kill w%maxstep%
|
149 |
|
|
kill wsq%maxstep%
|
150 |
|
|
kill n%maxstep%
|
151 |
|
|
kill vg%maxstep%
|
152 |
|
|
|
153 |
|
|
/* smooth, refine and combine each layer in turn
|
154 |
|
|
|
155 |
|
|
|
156 |
|
|
copy h%maxstep% hs%maxstep%
|
157 |
|
|
copy v%maxstep% vs%maxstep%
|
158 |
|
|
kill h%maxstep%
|
159 |
|
|
kill v%maxstep%
|
160 |
|
|
setcell hs%maxstep%
|
161 |
|
|
setwindow hs%maxstep%
|
162 |
|
|
|
163 |
|
|
&do j := %maxstep% &to 1 &by -1
|
164 |
|
|
&sv i [calc %j% - 1]
|
165 |
|
|
|
166 |
|
|
&type Refine from %j% to %i%
|
167 |
|
|
|
168 |
|
|
/* for the first stage where the coarser grid is refined and smoothed, set window to the coarse grid
|
169 |
|
|
setcell h%i%
|
170 |
|
|
setwindow maxof
|
171 |
|
|
|
172 |
|
|
/* create smoothed higher resolution versions of h and v_bg, hopefully with no nulls!
|
173 |
|
|
hs%j%_%i% = focalmean(hs%j%, circle, 2)
|
174 |
|
|
vs%j%_%i% = focalmean(vs%j%, circle, 2)
|
175 |
|
|
|
176 |
|
|
setcell h%i%
|
177 |
|
|
&describe h%i%
|
178 |
|
|
&sv cellsize %grd$dx%
|
179 |
|
|
&describe bboxgrid
|
180 |
|
|
setwindow [calc %grd$xmin% - 4 * %cellsize%] [calc %grd$ymin% - 4 * %cellsize%] [calc %grd$xmax% + 4 * %cellsize%] [calc %grd$ymax% + 4 * %cellsize%] h%i%
|
181 |
|
|
|
182 |
|
|
/* create no-null version of finer h and v
|
183 |
|
|
h%i%_c = con(isnull(h%i%), 0, h%i%)
|
184 |
|
|
v%i%_c = con(isnull(v%i%), %bigvar%, v%i%)
|
185 |
|
|
|
186 |
|
|
/* combine two values using least variance
|
187 |
|
|
hs%i% = (h%i%_c / v%i%_c + hs%j%_%i% / vs%j%_%i% ) / (1.0 / v%i%_c + 1.0 / vs%j%_%i%)
|
188 |
|
|
vs%i% = 1 / (1.0 / v%i%_c + 1.0 / vs%j%_%i%)
|
189 |
|
|
|
190 |
|
|
kill v%i%_c
|
191 |
|
|
kill h%i%_c
|
192 |
|
|
kill v%i%
|
193 |
|
|
kill h%i%
|
194 |
|
|
kill vs%j%_%i%
|
195 |
|
|
kill hs%j%_%i%
|
196 |
|
|
kill hs%j%
|
197 |
|
|
kill vs%j%
|
198 |
|
|
&end
|
199 |
|
|
|
200 |
|
|
/* result is hs0, with variance vs0
|
201 |
|
|
|
202 |
|
|
kill bboxgrid
|