1
|
&args DEMgrid
|
2
|
|
3
|
/* Calculates a measure of the noisiness of the surface.
|
4
|
/* John Gallant 8 Jan 2009
|
5
|
/*
|
6
|
/* Modified from noisemag on 15 Oct 2010 to work at 3 second resolution
|
7
|
/*
|
8
|
/* The principle is that noise is characterised by short-range variation.
|
9
|
/* The difference from the local mean will vary erratically over a short
|
10
|
/* distance in the presence of noise. Terrain features are generally smoother
|
11
|
/* so do not vary as erratically. The noise is measured as the local standard
|
12
|
/* deviation of the difference from the mean.
|
13
|
/*
|
14
|
/* The local mean is calculated using an annulus so that local peaks and
|
15
|
/* holes are further accentuated.
|
16
|
/*
|
17
|
/* The raw noise magnitude is smoothed using two steps of median filtering
|
18
|
/* to give the median value over a 25-cell circle.
|
19
|
/*
|
20
|
/* The calculated values can be used as the standard deviation of the
|
21
|
/* noise of the surface.
|
22
|
/*
|
23
|
/* One flaw is that higher relief terrain is considered "noisy" by this method.
|
24
|
/* (Lower relief terrain is not.)
|
25
|
/*
|
26
|
/* Could remove the effect of steeper terrain by using a fitted polynomial
|
27
|
/* rather than a simple mean as the basis for calculating the difference.
|
28
|
/* Not pursued at this time because the application is to control where
|
29
|
/* vegetation edges are visible in SRTM data, and high relief areas also
|
30
|
/* obscure the veg edges, so including the high relief as "high noise"
|
31
|
/* is not a problem.
|
32
|
|
33
|
setwindow %DEMgrid%
|
34
|
setcell %DEMgrid%
|
35
|
|
36
|
noisemean = focalmean(%DEMgrid% , annulus, 1, 2)
|
37
|
noisediffmn = %DEMgrid% - noisemean
|
38
|
noiseraw = focalstd(noisediffmn, circle, 2)
|
39
|
/* setcell minof
|
40
|
/* noisemag2 = focalmedian(aggregate(noiseraw, 2, median), circle, 2)
|
41
|
/* setcell %DEMgrid%
|
42
|
/* noisemag = resample(noisemag2, #, bilinear)
|
43
|
|
44
|
/* testing rejection of terrain signal by capturing longer-range differences from mean
|
45
|
/* noisemag_d3 looks particularly good - terrain effects almost completely absent, for both small and extended features
|
46
|
|
47
|
noisemndf = focalmean(noisediffmn, circle, 2) /* calculate mean difference from mean
|
48
|
noisemndf_abs = sqrt(1 + sqr(noisemndf)) - 1 /* soft absolute value of that
|
49
|
noisemndf_mdn = focalmedian(noisemndf_abs, circle, 3) /* smooth, this is the terrain signal
|
50
|
noiseraw_rdc = noiseraw / (1 + 2 * mndiffmn2bfmd) /* reduce noiseraw by the terrain signal
|
51
|
|
52
|
setcell minof
|
53
|
noisemag_2_2 = focalmedian(aggregate(noiseraw_rdc, 2, median), circle, 2)
|
54
|
setcell %DEMgrid%
|
55
|
noisemag = resample(noisemag_2_2, #, bilinear)
|