Project

General

Profile

Download (36.6 KB) Statistics
| Branch: | Revision:
1 3c841968 Jim Regetz
% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-
2
%
3
% Using Gregor Gorjanc's bash script:
4
%   sweave -old=<pdfviewer> <filename.Rnw>
5
%
6
% Using my simple bash script wrapper:
7
%   sweave_jr <basename>
8
%
9
% Using the lower level commands, in two steps: 
10
%   R CMD Sweave <filename.Rnw>
11
%   pdflatex <filename.tex>
12
%
13
% To extract R code from within R:
14
%   Stangle("<filename.Rnw>", annotate=FALSE)
15 e6e6f64f Jim Regetz
\documentclass{article}
16 8ceb6da6 Jim Regetz
\usepackage[font=small,labelfont=bf]{caption}
17 3c841968 Jim Regetz
\usepackage[font=normalsize,singlelinecheck=false]{subfig}
18 e6e6f64f Jim Regetz
\usepackage{enumitem}
19 3c841968 Jim Regetz
\usepackage[utf8]{inputenc}
20
\usepackage{a4wide}
21
\usepackage{verbatim}
22 0ffc32ec Jim Regetz
\usepackage{listings}
23 3c841968 Jim Regetz
\usepackage[colorlinks=true,urlcolor=red]{hyperref}
24 0ffc32ec Jim Regetz
\usepackage{color}
25 8ceb6da6 Jim Regetz
\usepackage{xspace}
26 3c841968 Jim Regetz
27 e6e6f64f Jim Regetz
% shrink list item spacing
28
\setlist{noitemsep}
29 8ceb6da6 Jim Regetz
30 0ffc32ec Jim Regetz
% set some listings options
31
\lstset{
32 e6e6f64f Jim Regetz
  basicstyle=\footnotesize\ttfamily,
33 0ffc32ec Jim Regetz
  stringstyle=\ttfamily,
34
  commentstyle=\itshape\ttfamily,
35
  showstringspaces=false,
36
}
37
38 e6e6f64f Jim Regetz
% define degree command for convenience
39
\newcommand{\N}{\ensuremath{^\circ\mbox{N}}\xspace}
40
\newcommand{\W}{\ensuremath{^\circ\mbox{W}}\xspace}
41
42
% add page numbers, but no headers
43
\thispagestyle{plain}
44
45
\title{SRTM/ASTER boundary analysis}
46 d126363c Jim Regetz
\author{Jim Regetz, NCEAS}
47
\date{Last update: 13 Jul 2011}
48 e6e6f64f Jim Regetz
49 3c841968 Jim Regetz
\begin{document}
50
51 8ceb6da6 Jim Regetz
\maketitle
52
53
% don't number document sections
54
\setcounter{secnumdepth}{-1}. 
55
56 d126363c Jim Regetz
% first compute some values to use in the text
57
<<echo=FALSE,results=hide>>=
58
  # Gaussian weighted average parameters
59
  gwa.r <- 0.001
60
  gwa.50.pix <- round(sqrt(-log(0.5)/gwa.r))
61
  gwa.50.km <- gwa.50.pix * 90 / 1000
62
  gwa.01.pix <- round(sqrt(-log(0.01)/gwa.r))
63
  gwa.01.km <- gwa.01.pix * 90 / 1000
64
65
  # mean/median ASTER-SRTM elevation differences
66
  library(raster)
67
  demdir <- "/home/regetz/media/temp/terrain/dem/"
68
  d.srtm <- raster(file.path(demdir, "srtm_150below.tif"))
69
  d.aster <- raster(file.path(demdir, "aster_300straddle.tif"))
70
  d.delta.vals <- values(crop(d.aster, extent(d.srtm))) - values(d.srtm)
71
  delta.median <- median(d.delta.vals)
72
  delta.mean <- mean(d.delta.vals)
73
  delta.sd <- sd(d.delta.vals)
74
  delta.q <- quantile(d.delta.vals, c(0.1, 0.25, 0.75, 0.9))
75
@
76
77
\paragraph{Brief summary of findings}
78 8ceb6da6 Jim Regetz
\begin{itemize}
79 d126363c Jim Regetz
 \item SRTM vs ASTER differences
80 8ceb6da6 Jim Regetz
  \begin{itemize}
81 d126363c Jim Regetz
   \item ASTER is systematically lower, by 12 meters in the median case
82
   \item \ldots but with variability: standard deviation of
83
     $\mbox{ASTER}_i-\mbox{SRTM}_i$ is \Sexpr{round(delta.sd, 1)}
84
   \item ASTER also has numerous spurious spikes
85 8ceb6da6 Jim Regetz
   \item ASTER has more high-frequency variability (``texture''),
86
     affecting slope/aspect?
87
  \end{itemize}
88 d126363c Jim Regetz
 \item Fusion via northward exponential rampdown of boundary delta
89 8ceb6da6 Jim Regetz
  \begin{itemize}
90 d126363c Jim Regetz
   \item eliminates elevation cliff at 60\N
91 8ceb6da6 Jim Regetz
   \item leaves abrupt transition in SRTM/ASTER textural differences
92 d126363c Jim Regetz
   \item introduces north-south ridging artifacts
93
   \item (\emph{no further treatment in this document})
94 8ceb6da6 Jim Regetz
  \end{itemize}
95 d126363c Jim Regetz
 \item Fusion via multiresolution spline
96 8ceb6da6 Jim Regetz
  \begin{itemize}
97 d126363c Jim Regetz
   \item eliminates elevation cliff at 60\N
98
   \item leaves abrupt transition in derived slope and aspect
99
   \item unclear whether derived values of aspect and flow direction
100
     in the transition zone are acceptable
101 8ceb6da6 Jim Regetz
  \end{itemize}
102 d126363c Jim Regetz
 \item Fusion via Gaussian weighted average of SRTM/ASTER
103 8ceb6da6 Jim Regetz
  \begin{itemize}
104 d126363c Jim Regetz
   \item eliminates elevation cliff at 60\N
105
   \item also smooths transition slope and aspect
106
   \item unclear whether derived values of aspect and flow direction in
107
     the transition zone are acceptable
108 8ceb6da6 Jim Regetz
  \end{itemize}
109
 \item Canada DEM itself has problems
110
  \begin{itemize}
111
   \item 60\N coincides with provincial boundaries; there are
112
     clear 60\N artifacts in this layer!
113
   \item other evident tiling artifacts too
114
  \end{itemize}
115
 \item Other comments
116
  \begin{itemize}
117
   \item N/S bias to aspect, flow direction computed on unprojected data
118
     at higher latitudes?
119
  \end{itemize}
120
\end{itemize}
121
122
\paragraph{To do (possibly)}
123
\begin{itemize}
124 d126363c Jim Regetz
 \item add constant offset of \Sexpr{-delta.median}m to ASTER
125
 \item apply low pass filter to ASTER to reduce high frequency
126
 variation?
127
 \item apply algorithm to remove spikes (\ldots but maybe beyond scope?)
128 8ceb6da6 Jim Regetz
\end{itemize}
129 3c841968 Jim Regetz
130 8ceb6da6 Jim Regetz
\clearpage
131 0ffc32ec Jim Regetz
132 e6e6f64f Jim Regetz
%-----------------------------------------------------------------------
133 8ceb6da6 Jim Regetz
\section{Terrain layer production methodology}
134 e6e6f64f Jim Regetz
%-----------------------------------------------------------------------
135 8ceb6da6 Jim Regetz
136 d126363c Jim Regetz
For the purposes of assessing artifacts associated with the northern
137
boundary between SRTM and ASTER, I focused on a narrow band along the
138
60\N boundary in Canada (Figure \ref{focal-area}). The latitudinal
139
extent of this band is 59.875\N - 60.125\N (i.e., 300 3" pixels
140
straddling 60\N), and the longitudinal extent is 136\W to 96\W (i.e.,
141
48,000 pixels wide). See Listing \ref{code-gdalwarp} for code. Within
142
this focal region, I generated latitudinal profiles of mean elevation,
143
slope, aspect, and flow direction, using the separate SRTM and ASTER
144
component DEMs themselves, using several different fused ASTER/SRTM DEMs
145
(see below), and using the Canadian Digital Elevation Data (CDED) as an
146
independent reference layer
147
(\url{http://www.geobase.ca/geobase/en/data/cded/index.html}). I also
148
then computed latitudinal correlations and RMSEs between the fused
149
layers and each of SRTM, ASTER, and CDED.
150 8ceb6da6 Jim Regetz
151
\paragraph{Elevation} This document includes latitudinal
152
characterizations of terrain values based on three different variants of
153 d126363c Jim Regetz
a fused 3" ASTER/SRTM DEM. I also explored (and briefly describe below) two
154 8ceb6da6 Jim Regetz
additional approaches to fusing the layers, but do not include further
155
assessment of these here.
156 0ffc32ec Jim Regetz
157 9fd6a7de Jim Regetz
\subparagraph{Simple fusion} Naive concatenation of SRTM below 60\N with
158 8ceb6da6 Jim Regetz
ASTER above 60\N, without applying any modifications to deal with
159
boundary artifacts (Figure \ref{blend-simple}).
160 0ffc32ec Jim Regetz
161 8ceb6da6 Jim Regetz
\subparagraph{Multiresolution spline} Application of Burt \& Adelson's
162
(1983) method for blending overlapping images using multiresolution
163 d126363c Jim Regetz
splines, as implemented in the \emph{Enblend} software package (version
164
4.0, \url{http://enblend.sourceforge.net}). Data preparation and
165
post-processing were handled in R (see Listing \ref{code-enblend}). As
166
presented here, the SRTM and ASTER inputs were prepared such that the
167
overlap zone is 75 latitudinal rows (6.75km) (Figure
168
\ref{blend-multires}).
169 0ffc32ec Jim Regetz
170 8ceb6da6 Jim Regetz
\subparagraph{Gaussian weighted average} Blend of the two layers using
171
weighted averaging such that the relative contribution of the SRTM
172 d126363c Jim Regetz
elevation is zero at 60\N, and increases as a function of distance
173
moving south away from 60\N (Equation \ref{eq-gaussian}).
174 8ceb6da6 Jim Regetz
\begin{eqnarray}
175
\label{eq-gaussian}
176
&fused_{x,y} =
177
  \left\{
178
    \begin{array}{c l}
179
      ASTER_{x,y} & \mbox{above } 60^\circ \mbox{N} \\
180
      w_{x,y}ASTER_{x,y} + (1-w_{x,y})SRTM_{x,y} & \mbox{below } 60^\circ \mbox{N}
181
    \end{array}
182
  \right.\\
183
&\mbox{where }
184
  w_{x,y}=e^{-rD_{y}^{2}}
185
  \mbox{ and }
186
  D_{y} \mbox{ is the distance from } 60^\circ \mbox{N in units of pixels.} \nonumber
187
\end{eqnarray}
188 d126363c Jim Regetz
For the assessment presented here, the weighting function function was
189
parameterized using $r$=\Sexpr{gwa.r}, producing equal weights for SRTM
190
and ASTER at a distance of $\sim$\Sexpr{round(gwa.50.km, 1)}km
191
(\Sexpr{gwa.50.pix} cells) south of the the boundary, and a relative
192
weight for ASTER of only 1\% by $\sim$\Sexpr{round(gwa.01.km, 1)}km
193
(\Sexpr{gwa.01.pix} cells) (Figure \ref{blend-gaussian}). See
194
\texttt{OPTION 3} in Listing \ref{code-correct}.
195
196
\subparagraph{Others not shown} I also experimented with some additional
197
fusion approaches, but have excluded them from further analysis in this
198
document.
199
200
\emph{Fused with exponential ramp north of 60\N.}
201
The first step was to take the pixel-wise difference between SRTM and
202 9fd6a7de Jim Regetz
ASTER in the row immediately below 60\N (i.e., the northernmost extent of
203 d126363c Jim Regetz
SRTM). An exponentially declining fraction of this difference was then
204
then added back into the ASTER values north of 60\N. This does a fine
205
job of eliminating the artificial shelf and thus the appearance of a
206
seam right along the 60\N boundary, but it does not address the abrupt
207 8ceb6da6 Jim Regetz
transition in texture (i.e., the sudden appearance of high frequency
208
variability moving north across the boundary). Additionally, it
209
introduces vertical ``ridges'' running north from the boundary. These
210
arise because the calculated ramps are independent from one longitudinal
211
``column'' to the next, and thus any changes in the boundary difference
212
from one pixel to the next lead to adjacent ramps with different
213
inclines.
214 0ffc32ec Jim Regetz
215 d126363c Jim Regetz
\emph{Simple LOESS predictive model.}
216
This involved first calculating the difference between SRTM and ASTER
217
everywhere south of 60\N, and then fitting a LOESS curve to these
218
differences using the actual ASTER elevation as a predictor. I then used
219
the fitted model to predict the ASTER-SRTM difference for each ASTER
220
cell north of 60\N, and added a declining fraction (based on a Gaussian
221
curve) of this difference to the corresponding ASTER elevation.
222
Conceptually, this amounts to applying an ASTER-predicted SRTM
223
correction to the ASTER elevation values, where the correction term has
224
a weight that declines to zero with increasing distance (north) away
225
from the boundary. However, this method alone didn't yield particularly
226
promising results in removing the 60\N seam itself, presumably because
227
adding a predicted correction at the boundary does not close the
228
SRTM-ASTER gap nearly as efficiently as do corrections based directly on
229
the observed SRTM vs ASTER elevation differences. I therefore haven't
230 8ceb6da6 Jim Regetz
pursued this any further, although it (or something like it) may prove
231
useful in combination with one of the other methods.
232 0ffc32ec Jim Regetz
233 8ceb6da6 Jim Regetz
\paragraph{Slope}
234
For each of the three main fused DEM variants described above, slope was
235 d126363c Jim Regetz
calculated using \texttt{gdaldem} (GDAL 1.8.0, released 2011/01/12):
236 8ceb6da6 Jim Regetz
\begin{verbatim}
237
    $ gdaldem slope -s 111120 <input_elevation> <output_slope>
238
\end{verbatim}
239
Note that the scale option used here is as recommended in the
240
\texttt{gdaldem} documentation:
241
\begin{quote}
242 d126363c Jim Regetz
  ``\emph{If the horizontal unit of the source DEM is degrees (e.g
243
  Lat/Long WGS84 projection), you can use scale=111120 if the vertical'
244
  units are meters}''
245 8ceb6da6 Jim Regetz
\end{quote}
246
The output slope raster is in units of degrees.
247
248
\paragraph{Aspect}
249
As was the case with slope, aspect was calculated using
250
\texttt{gdaldem}:
251
\begin{verbatim}
252
    $ gdaldem aspect -s 111120 <input_elevation> <output_aspect>
253
\end{verbatim}
254
The output aspect raster values indicate angular direction in units of
255
degrees, with 0=North and proceeding clockwise.
256
257
\paragraph{Flow direction}
258 d126363c Jim Regetz
259
Flow direction was calculated using the GRASS (GRASS GIS 6.4.1)
260
\texttt{r.terraflow} module; see code listing \ref{code-flowdir}.
261
Because of a $\sim$30k ($2^{15}$) limit to the input raster dimension
262
size in the pre-built GRASS \texttt{r.terraflow} module I used, this
263
analysis was restricted to a smaller longitudinal subset of the data,
264
spanning 125\W to 100\W.
265
266
The default flow direction output of this module is encoded so as to
267
indicate \emph{all} downslope neighbors, also known as the Multiple Flow
268
Direction (MFD) model. However, to simplify post-processing and
269
summarization, the results here are based on an alternative Single Flow
270
Direction (SFD, \emph{a.k.a.} D8) model, which indicates the neighbor
271
associated with the steepest downslope gradient. Note that this is
272
equivalent to what ArcGIS GRID \texttt{flowaccumulation} command does. I
273
then recoded the output raster to use the same azimuth directions used
274
by \texttt{gdaldem aspect}, as described above for aspect.
275 0ffc32ec Jim Regetz
276
277 e6e6f64f Jim Regetz
%-----------------------------------------------------------------------
278
\section{Latitudinal mean terrain profiles}
279
%-----------------------------------------------------------------------
280
281 8ceb6da6 Jim Regetz
282 d126363c Jim Regetz
\paragraph{Elevation} SRTM, ASTER, and CDED all share a very
283 8ceb6da6 Jim Regetz
similarly shaped mean elevation profile, but with differing heights
284
(Figure \ref{mean-elevation}). SRTM tends to be highest, ASTER is
285 d126363c Jim Regetz
lowest, and CDED is intermediate. The magnitude of average difference
286
between SRTM and ASTER is fairly consistent not only across latitudes,
287
but also across elevations (Figure \ref{aster-srtm}). The overall median
288
difference between ASTER and SRTM (i.e., considering
289
$\mbox{ASTER}_i-\mbox{SRTM}_i$ for all pixels $i$ where the two DEMs
290
co-occur) is \Sexpr{delta.median} meters, with a mean of
291 8711607d Jim Regetz
\Sexpr{round(delta.mean, 2)} meters, and this more or less holds (within
292
a few meters) across the observed range of elevations (Figure
293 d126363c Jim Regetz
\ref{aster-srtm-bins}). However, while this average offset is broadly
294
consistent across latitudes and across elevation zones, additional
295
variation is evident at the pixel level. Again focusing on the
296
pixel-wise differences, they appear to be symmetrically distributed
297
about the mean with a standard deviation of \Sexpr{round(delta.sd, 1)}
298
and quartiles ranging from \Sexpr{delta.q["25%"]} to
299
\Sexpr{delta.q["75%"]} meters; ASTER elevations are actually greater
300
than SRTM for \Sexpr{round(100*mean(d.delta.vals>0))}\% of pixels (see
301
Figure \ref{aster-srtm-scatter}). Thus, although adding a constant
302
offset of \Sexpr{-delta.median} meters to the ASTER DEM would clearly
303 9fd6a7de Jim Regetz
center it with respect to the SRTM (at least in the Canada focal
304 d126363c Jim Regetz
region), appreciable differences would remain. Figure
305
\ref{aster-srtm-scatter} also highlights the existence of several
306
obviously spurious ASTER spikes of >1000m; although not shown here,
307
these tend to occur in small clumps of pixels, perhaps corresponding to
308
false elevation readings associated with clouds?
309 8ceb6da6 Jim Regetz
310
Not surprisingly, simple fusion produces an artificial $\sim$12m cliff
311
in the mean elevation profile (Figure \ref{mean-elevation}). At least in
312
terms of mean elevation, this artifact is completely removed by both the
313 9fd6a7de Jim Regetz
multiresolution spline and Gaussian weighted average methods. The
314 8ceb6da6 Jim Regetz
transition is, to the eye, slightly smoother in the former case,
315
although ultimately this would depend on the chosen zone of overlap and
316
on the exact parameterization of the weighting function.
317
318
\paragraph{Slope} The mean ASTER slope is uniformly steeper than the
319
mean SRTM slope at all latitudes in the area of overlap, by nearly 1
320
degree (Figure \ref{mean-slope}). However, the shape of the profile
321
itself is nearly identical between the two. Although this may partly
322
reflect inherent SRTM vs ASTER differences, my guess is that CGIAR
323 9fd6a7de Jim Regetz
post-processing of the particular SRTM product we're using has removed
324 8ceb6da6 Jim Regetz
some of the high frequency ``noise'' that remains in ASTER?
325
\textbf{\color{red}[todo: check!]}
326
327 d126363c Jim Regetz
Note that the CDED tends to be flatter than both SRTM and
328
ASTER (presumably because it is at least partially derived from
329
contour-based data \textbf{\color{red}[todo: check!]}). Moreover, this
330
figure makes it clear that CDED has some major artifacts at
331 e6e6f64f Jim Regetz
regular intervals. The spike especially at 60\N (which coincides with
332
provincial boundaries across the entirety of western Canada) means we
333 d126363c Jim Regetz
probably need to scuttle our plans to use this DEM as a formal reference
334 e6e6f64f Jim Regetz
dataset for boundary analysis.
335 0ffc32ec Jim Regetz
336
The simple fused layer exhibits a dramatic spike in slope at the
337 8ceb6da6 Jim Regetz
immediate 60\N boundary, undoubtedly associated with the artificial
338 0ffc32ec Jim Regetz
elevation cliff. This artifact is eliminated by both the multiresolution
339 d126363c Jim Regetz
spline and Gaussian weighting. However, the former exhibits a sudden step
340 0ffc32ec Jim Regetz
change in slope in the SRTM-ASTER overlap region, whereas the transition
341
is smoothed out in the latter. This likely reflects the fact that the
342
multiresolution spline effectively uses a very narrow transition zone
343
for stitching together high frequency components of the input images,
344 8ceb6da6 Jim Regetz
and it seems likely that these are precisely the features responsible
345
for the shift in mean slope. 
346
347
\paragraph{Aspect}
348
For the purposes of latitudinal profiles, aspect values were summarized
349
using a circular mean (Equation \ref{eq-cmean}).
350
\begin{equation}
351
 \bar{x} = \mathrm{atan2}
352
   \left(
353
    \sum_{i=1}^{n}\frac{\sin(x_i)}{n},
354
    \sum_{i=1}^{n}\frac{\cos(x_i)}{n}
355
   \right )
356
\label{eq-cmean}
357
\end{equation}
358
where $x_i$ is the aspect value (in radians) of pixel $i$.
359
360
As indicated in Figure \ref{mean-aspect}, the circular mean aspect
361 d126363c Jim Regetz
values of SRTM and ASTER are generally similar across all latitudes in
362
the area of overlap, and mean aspect values calculated on CDED
363 8ceb6da6 Jim Regetz
are similar at most but not all latitudes. The mean values at nearly all
364
latitudes are directed either nearly north or nearly south, though
365
almost always with a slight eastward rather than westward inclination.
366
In fact, there appears to be a general bias towards aspect values
367
orienting along the north-south axis, as is especially apparent in the
368
rose diagrams of Figure \ref{rose-diag-aspect}. I suspect this is an
369
artifact of our use of unprojected data, especially at these high
370
latitudes. Because the unprojected raster is effectively 'stretched'
371
east and west relative to the actual topography, elevational gradients
372
along the east-west axis are artificially flattened, and the direction
373
of dominant gradient is more likely to be along the north-south axis.
374
375
Upon further reflection, it's not clear whether these patterns of mean
376
aspect are particular useful diagnostics, as they seems to be sensitive
377
to subtle variations in the data. Referring again to Figure
378
\ref{rose-diag-aspect}, note how the mean direction flips from nearly
379
north to nearly south between the two latitudes, even though the
380
distributions of pixel-wise aspect values are nearly indistinguishable
381
by eye.
382
383
In any case, not surprisingly, the simple fused layer matches the SRTM
384
aspect values south of 60\N and the ASTER aspect values north of 60\N; at
385
the immediate boundary, the mean aspect is northward, as one would
386
expect in the presence of a cliff artifact at the seam.
387
388
Interestingly, the aspect layers derived from the two blended DEMs
389 9fd6a7de Jim Regetz
(multiresolution spline and Gaussian weighted average) exhibit a
390 8ceb6da6 Jim Regetz
consistent mean northward inclination at all latitudes in their
391
respective fusion zones. This pattern is visually obvious at latitudes
392
between 59.95\N and 60\N in the bottom two panels of Figure
393
\ref{mean-aspect}. This is almost certainly a signal of the blending of
394
the lower elevation ASTER to the north with higher elevation SRTM to the
395 d126363c Jim Regetz
south, introducing a north-facing tilt (however slight) to the data
396
throughout this zone. 
397 8ceb6da6 Jim Regetz
398 e6e6f64f Jim Regetz
\paragraph{Flow direction} With the exception of edge effects at the
399 d126363c Jim Regetz
margins of the input rasters, mean ASTER-derived flow is nearly
400
northward at all latitudes, and SRTM-derived flow is nearly northward at
401
all but a few latitudes (Figure \ref{mean-flowdir}). This seems
402
reasonable considering that most pixels in this Canada test region fall
403
in the Arctic drainage. For unknown reasons, CDED produces southward
404
mean flow direction at numerous latitudes, and generally seems to have a
405
slightly more eastward tendency. As was the case with aspect, note that
406
a general north-south bias is evident (Figure \ref{rose-diag-flowdir}),
407
again likely due to use of an unprojected raster at high latitudes.
408 e6e6f64f Jim Regetz
409
The various fused layer profiles look as one would expect (Figure
410
\ref{mean-flowdir}), although the overall lack of latitudinal
411
variability in mean flow direction in SRTM, ASTER, and all three derived
412
layers makes it hard to say much more than that. 
413
414
%-----------------------------------------------------------------------
415
\section{Informal correlation analysis}
416
%-----------------------------------------------------------------------
417
418
\paragraph{Elevation}
419 d126363c Jim Regetz
Pearson correlations between SRTM and ASTER are quite high, typically
420
>0.999, and RMSEs are on the order of 10-15 meters (Figures
421
\ref{corr-elevation} and \ref{rmse-elevation}). Spikes (downward for
422
correlation, upward for RMSE) occur at some latitudes, quite possibly
423
associated with the observed extreme spikes in the ASTER DEM itself. As
424
expected, the multiresolution spline and Gaussian weighted average both
425
produce layers that gradually become less similar to ASTER and more
426
similar to SRTM moving south from 60\N, but in slightly different ways.
427
This gradual transition is less evident when considering associations
428
with CDED (bottom panels of Figures \ref{corr-elevation} and
429 e6e6f64f Jim Regetz
\ref{rmse-elevation}), which in general is less correlated with SRTM,
430 d126363c Jim Regetz
and even less with ASTER, than those two layers are with each other.
431 e6e6f64f Jim Regetz
432
\paragraph{Slope}
433
The patterns of slope layer similarity are much like those described
434
above for elevation, although the correlations are somewhat lower
435
($\sim$0.94 between SRTM and ASTER (Figure \ref{corr-slope})). RMSEs
436
between SRTM and ASTER are approximately 2 at all latitudes where the
437
data co-occur (Figure \ref{rmse-slope}). Another difference, echoing a
438
pattern previously noted in the profiles of mean slope itself, is that
439
Gaussian weighted averaging produces a layer that exhibits a gradual
440
transition from ASTER to SRTM, whereas the multiresolution spline yields
441
an abrupt transition. Not surprisingly, the simple fused layer is even
442 9fd6a7de Jim Regetz
worse, producing not only a sudden transition but also aberrant values
443 d126363c Jim Regetz
at the fusion seam itself; note downward (upward) spikes in correlation
444 e6e6f64f Jim Regetz
(RMSE) at 60\N in the first column of plots in Figures \ref{corr-slope}
445
and \ref{rmse-slope}.
446
447
\paragraph{Aspect}
448
Because aspect values are on a circular scale, I calculated modified
449
versions of the Pearson correlation coefficient using the
450
\textbf{circular} R package function \texttt{cor.circular}. I believe
451
this implements the formula described by Jammalamadaka \& Sarma (1988):
452
\begin{equation}
453
r = \frac
454
  {\sum_{i=1}^{n} \sin(x_i-\bar{x}) \sin(y_i-\bar{y})}
455
  {\sqrt{\sum_{i=1}^{n} \sin(x_i-\bar{x})^2}
456
   \sqrt{\sum_{i=1}^{n} \sin(y_i-\bar{y})^2}}
457
\label{eq-ccor}
458
\end{equation}
459
where $x_i$ and $y_i$ are aspect values (in radians) for pixel $i$, and
460
$\bar{x}$ and $\bar{y}$ are circular means calculated as in Equation
461
\ref{eq-cmean}.
462
463 d126363c Jim Regetz
To calculate RMSEs for aspect, I did not attempt to use trigonometric
464
properties analogous to computation of the circular correlation, but
465
instead just imposed a simple correction whereby all pairwise
466
differences were computed using the shorter of the two paths around the
467
compass wheel (Equation \ref{eq-circ-rmse}). For example, the difference
468
between 0$^\circ$ and 150$^\circ$ is 150$^\circ$, but the difference
469
between 0$^\circ$ and 250$^\circ$ is 110$^\circ$.
470 e6e6f64f Jim Regetz
471
\begin{eqnarray}
472
\label{eq-circ-rmse}
473
&RMSE = \sqrt{\frac{\sum_{i=1}^{n} (\Delta_i^2)}{n}}\\ \nonumber
474
&\mbox{where }
475
  \Delta_i = \mbox{argmin} \left (|x_i-y_i|, 360-|x_i - y_i| \right )
476
\end{eqnarray}
477
478
Circular correlations between SRTM and ASTER were surprisingly low,
479
typically hovering around 0.5, but dipping down towards zero at numerous
480
latitudes (Figure \ref{corr-flowdir}). The corresponding RMSE is close
481
to 70 at all latitudes, surprisingly high considering that the maximum
482
difference between any two pixels is 180$^\circ$ (Figure
483 d126363c Jim Regetz
\ref{rmse-flowdir}). Comparison of SRTM and ASTER with CDED
484 e6e6f64f Jim Regetz
yields similar patterns.
485
486
For both of the blended layers (multiresolution spline and Gaussian
487
weighted average), aspect values calculated in the zone of ASTER/SRTM
488
overlap appear to be \emph{less} similar to the component DEMs than
489
those to data sources are to each other. As evident in the upper middle
490
and upper right panels of Figure \ref{corr-flowdir}, correlations
491
between fused aspect and aspect based on the original SRTM and ASTER
492
images are substantially \emph{negative} for many latitudes in the zone
493
of overlap. I don't have a good feel for what's going on here, although
494 d126363c Jim Regetz
it may just involve properties of the circular correlation statistic
495
that I don't have a good feel for. Note that the latitudinal profile of
496
aspect RMSEs is less worrisome (upper middle and upper right panels of
497
Figure \ref{rmse-flowdir}) and more closely resembles the profile of
498
slope RMSEs discussed above, particularly as regards the more gradual
499 e6e6f64f Jim Regetz
transition in case of Gaussian weighted averaging compared to the
500
multiresolution spline.
501
502
\paragraph{Flow direction}
503
As with aspect, circular correlation coefficients and adjusted RMSEs
504 d126363c Jim Regetz
were calculated for each latitudinal band. Flow direction correlations
505 e6e6f64f Jim Regetz
between SRTM and ASTER are slightly lower than for aspect, typically
506
around 0.40 but spiking negative at a handful of latitudes (Figure
507 d126363c Jim Regetz
\ref{corr-flowdir}). RMSEs hover consistently around $\sim$75, slightly
508 e6e6f64f Jim Regetz
lower than was the case with aspect (Figure \ref{rmse-flowdir}).
509
510 d126363c Jim Regetz
Aside from an expected flow artifact at the southern image edge, and an
511 e6e6f64f Jim Regetz
unexpected (and as-yet unexplained) negative spike at $\sim$59.95\N, the
512
correlation profiles are fairly well-behaved for both blended layers,
513
and don't show the same odd behavior as was the case for aspect. Again
514 9fd6a7de Jim Regetz
it is clear that the multiresolution spline results in a much more abrupt
515 e6e6f64f Jim Regetz
transition than does the Gaussian weighted average. The RMSE flow
516
direction profiles echo this pattern (Figure \ref{rmse-flowdir}), and
517
indeed look almost the same as those computed using aspect (Figure
518
\ref{rmse-aspect}).
519
520
%-----------------------------------------------------------------------
521
% FIGURES
522
%-----------------------------------------------------------------------
523
524
\clearpage
525
526
%
527
% Figure: Boundary analysis region
528
%
529
\begin{figure}[h]
530 8711607d Jim Regetz
  \caption{Focal area used as the basis for boundary assessment. Note
531
  that for flow direction, analysis was restricted to a smaller
532
  longitudinal span (125\W to 100\W).}
533 e6e6f64f Jim Regetz
  \centering
534
<<echo=FALSE,results=hide,fig=TRUE,height=5,width=7>>=
535
  par(omi=c(0,0,0,0))
536
  library(raster)
537
  library(maps)
538
  demdir <- "/home/regetz/media/temp/terrain/dem/"
539
  dem <- raster(file.path(demdir, "fused_300straddle.tif"))
540
  map("world", xlim=c(-140,-70), ylim=c(40, 70), fill=TRUE, col="grey",
541
      mar=c(4,0,0,0))
542
  grid(col="darkgray")
543
  axis(1)
544
  axis(2)
545
  box()
546
  plot(extent(dem), col="red", add=TRUE)
547
  arrows(-110, 57, -113, 59.7, col="red", length=0.1)
548
  text(-110, 57, labels="focal region", col="red", font=3, pos=1)
549
@
550
  \label{focal-area}
551
\end{figure}
552
553
%
554
% Figure: Fusion methods
555
%
556
\begin{figure}
557
  \caption{SRTM/ASTER DEM fusion configurations}
558
  \centering
559
  \subfloat[Simple fusion]{
560
    \label{blend-simple}
561
<<echo=FALSE,results=hide,fig=TRUE, height=3, width=6>>=
562
  par(omi=c(0,0,0,0), mar=c(4,4,0,2)+0.1)
563
  plot(0, xlim=c(-150, 150), ylim=c(0, 1), type="n", xaxt="n", yaxt="n",
564
      bty="n", xlab="Latitude", ylab=NA)
565
  rect(0, 0, 150, 0.75, col="red", density=5, angle=45)
566
  rect(-150, 0, 0, 0.75, col="blue", density=5, angle=135)
567
#  axis(1, at=seq(-120, 120, by=60), labels=seq(59.9, 60.1, by=0.05))
568
  axis(1, at=seq(-150, 150, by=75), labels=c(59.875, NA, 60, NA, 60.125))
569
  text(-150, 0.9, labels="SRTM", pos=4, col="blue")
570
  text(150, 0.9, labels="ASTER", pos=2, col="red")
571
@
572
  }\\
573
  \subfloat[Multiresolution spline]{
574
    \label{blend-multires}
575
<<echo=FALSE,results=hide,fig=TRUE, height=3, width=6>>=
576
  par(omi=c(0,0,0,0), mar=c(4,4,0,2)+0.1)
577
  plot(0, xlim=c(-150, 150), ylim=c(0, 1), type="n", xaxt="n", yaxt="n",
578
      bty="n", xlab="Latitude", ylab=NA)
579
  rect(-75, 0, 0, 0.75, col="lightgray")
580
  rect(-75, 0, 150, 0.75, col="red", density=5, angle=45)
581
  rect(-150, 0, 0, 0.75, col="blue", density=5, angle=135)
582
  #axis(1, at=seq(-150, 150, by=60), labels=seq(59.875, 60.125, by=0.05))
583
  #axis(1, at=0, labels=60, col="red", cex=0.85)
584
  #axis(1, at=seq(-120, 120, by=60), labels=seq(59.9, 60.1, by=0.05))
585
  axis(1, at=seq(-150, 150, by=75), labels=c(59.875, NA, 60, NA, 60.125))
586
  text(-150, 0.9, labels="SRTM", pos=4, col="blue")
587
  text(-37.5, 0.8, labels="overlap", pos=3, font=3)
588
  text(150, 0.9, labels="ASTER", pos=2, col="red")
589
@
590
  }\\
591
  \subfloat[Gaussian weighted average]{
592
    \label{blend-gaussian}
593
<<echo=FALSE,results=hide,fig=TRUE, height=3, width=6>>=
594
  par(omi=c(0,0,0,0), mar=c(4,4,1,2)+0.1)
595
  curve(exp(-0.001*x^2), from=-150, to=0, xlim=c(-150, 150), col="red",
596
      xaxt="n", bty="n", xlab="Latitude", ylab="Weighting")
597
  curve(1+0*x, from=0, to=150, add=TRUE, col="red")
598
  curve(1-exp(-0.001*x^2), from=-150, to=0, add=TRUE, col="blue")
599
  #axis(1, at=seq(-150, 150, by=60), labels=seq(59.875, 60.125, by=0.05))
600
  #axis(1, at=0, labels=60, col="red", cex=0.85)
601
  #axis(1, at=seq(-120, 120, by=60), labels=seq(59.9, 60.1, by=0.05))
602
  axis(1, at=seq(-150, 150, by=75), labels=c(59.875, NA, 60, NA, 60.125))
603
  text(-100, 0.6, labels="gaussian\nblend")
604
  text(-140, 0.9, labels="SRTM", col="blue", pos=4)
605
  text(-140, 0.1, labels="ASTER", col="red", pos=4)
606
@
607
  }
608
\end{figure}
609
610
%
611
% Figure: Mean latitudinal profiles
612
%
613
\begin{figure}
614
  \centering
615
  \caption{Mean elevation (m)}
616
    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../dem/elevation-assessment.pdf}
617
    \label{mean-elevation}
618
\end{figure}
619
\begin{figure}
620
  \centering
621
  \caption{Mean slope (degrees)}
622
    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../slope/slope-assessment.pdf}
623
    \label{mean-slope}
624
\end{figure}
625
\begin{figure}
626
  \centering
627
  \caption{Circular mean aspect (azimuth)}
628
    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../aspect/aspect-assessment.pdf}
629
    \label{mean-aspect}
630
\end{figure}
631
\begin{figure}
632
  \centering
633
  \caption{Circular mean flow direction}
634
    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../flow/flowdir-assessment.pdf}
635
    \label{mean-flowdir}
636
\end{figure}
637
638
%
639
% Figure: ASTER vs SRTM elevation patterns
640
%
641
\begin{figure}
642
  \caption{ASTER vs SRTM elevation comparisons}
643
  \centering
644
  \subfloat[Boxplots of the arithmetic difference in elevations (ASTER -
645 d126363c Jim Regetz
  SRTM), summarized across a series of SRTM elevation bins. Boxes
646 e6e6f64f Jim Regetz
  include the median (horizontal band), 1st and 3rd quartiles (box
647 d126363c Jim Regetz
  extents), and $\pm$1.5$\times$IQR (whiskers); outliers not shown. Gray
648
  numbers above each whisker indicates how many thousands of pixels are
649
  included in the corresponding summary. Dashed red line indicates the
650
  median difference across all pixels south of 60\N.]{
651
    \includegraphics[width=\linewidth, trim=0 0 0 0.5in, clip=TRUE]{../dem/aster-srtm-bins.png}
652 e6e6f64f Jim Regetz
    \label{aster-srtm-bins}
653
  }\\
654 d126363c Jim Regetz
  \subfloat[Pixel-wise plot of SRTM vs ASTER for all values in the 150
655 e6e6f64f Jim Regetz
  latitudinal rows of overlap south of 60\N. Dashed blue line indicates
656
  the 1:1 diagonal, and the parallel red line is offset lower by the
657 8711607d Jim Regetz
  observed median difference between the two DEMs (\Sexpr{-delta.median}m).
658 d126363c Jim Regetz
  Inset histogram shows distribution of differences, excluding absolute
659 8711607d Jim Regetz
  differences >60m.]{
660 d126363c Jim Regetz
    \includegraphics[width=\linewidth, trim=0 0 0 0.5in, clip=TRUE]{../dem/aster-srtm-scatter.png}
661 e6e6f64f Jim Regetz
    \label{aster-srtm-scatter}
662
  }
663
  \label{aster-srtm}
664
\end{figure}
665
666
%
667
% Figures: Aspect rose diagrams
668
%
669 8ceb6da6 Jim Regetz
\begin{figure}[h!]
670
  \caption{Binned frequency distributions of aspect within two selected
671
    latitudinal strips: (a) within the SRTM-ASTER blend zone, and (b)
672
    south of the overlap zone. Aspect values here are based on the DEM
673
    produced via Gaussian weighted averaging, but similar patterns are
674
    evident using the other DEMs. Red spoke lines indicate the circular
675
    mean direction across all pixels at the associated latitude.}
676
  \centering
677
<<echo=FALSE,results=hide,fig=TRUE,height=5,width=7>>=
678
  library(circular)
679
  aspdir <- "/home/regetz/media/temp/terrain/aspect"
680
  a.bg <- raster(file.path(aspdir, "fused_300straddle_blendgau_a.tif"))
681
  par(mfrow=c(1,2), mar=c(0,0,4,0))
682
  y1 <- 200
683
  y2 <- 220
684
  cx <- circular(as.matrix(a.bg)[c(y1,y2),], units="degrees",
685
      rotation="clock", zero=pi/2)
686
  rose.diag(cx[1,], bins=8, axes=FALSE)
687
  mtext(paste("(a)", round(yFromRow(a.bg, y1), 3), "degrees"))
688
  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),
689
      labels=c("N", "W", "S", "E"))
690
  points(mean.circular(cx[1,], na.rm=TRUE), col="red")
691
  lines.circular(c(mean.circular(cx[1,], na.rm=TRUE), 0), c(0,-1),
692
      lwd=2, col="red")
693
  rose.diag(cx[2,], bins=8, axes=FALSE)
694
  mtext(paste("(b)", round(yFromRow(a.bg, y2), 3), "degrees"))
695
  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),
696
      labels=c("N", "W", "S", "E"))
697
  points(mean.circular(cx[2,], na.rm=TRUE), col="red")
698
  lines.circular(c(mean.circular(cx[2,], na.rm=TRUE), 0), c(0,-1),
699
      lwd=2, col="red")
700
@
701
  \label{rose-diag-aspect}
702
\end{figure}
703
704 e6e6f64f Jim Regetz
%
705
% Figures: Flow direction rose diagrams
706
%
707 8ceb6da6 Jim Regetz
\begin{figure}[h!]
708
  \caption{Relative frequencies of D8 flow directions (indicated by
709
    relative heights of the blue spokes) within two selected
710
    latitudinal bands: (a) within the SRTM-ASTER blend zone, and (b)
711
    south of the overlap zone. Flow directions here are based on the DEM
712
    produced via Gaussian weighted averaging. Red spoke lines indicate
713
    the circular mean flow direction across all pixels at the associated
714
    latitude.}
715
  \centering
716
<<echo=FALSE,results=hide,fig=TRUE,height=5,width=7>>=
717
  library(circular)
718
  flowdir <- "/home/regetz/media/temp/terrain/flow"
719
  # create function to recode terraflow SFD values into degrees, with
720
  # 0=North and proceeding clockwise (this matches gdaldem's default
721
  # azimuth output for aspect calculation)
722
  recode <- function(r) {
723
      v <- values(r)
724
      v[v==0] <- NA
725
      v[v==1] <- 90  ## east
726
      v[v==2] <- 135
727
      v[v==4] <- 180  ## south
728
      v[v==8] <- 225
729
      v[v==16] <- 270  ## west
730
      v[v==32] <- 315
731
      v[v==64] <- 0  ## north
732
      v[v==128] <- 45
733
      r[] <- v
734
      return(r)
735
  }
736
  sfd.bg <- recode(raster(file.path(flowdir,
737
      "fused_300straddle_blendgau_sfd.tif")))
738
  par(mfrow=c(1,2), mar=c(0,0,4,0))
739
  y1 <- 200
740
  y2 <- 220
741
  cx <- circular(as.matrix(sfd.bg)[c(y1,y2),], units="degrees",
742
      rotation="clock", zero=pi/2)
743
  rose.diag(cx[1,], bins=1000, border="blue", ticks=FALSE, axes=FALSE)
744
  mtext(paste("(a)", round(yFromRow(sfd.bg, y1), 3), "degrees"))
745
  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),
746
      labels=c("N", "W", "S", "E"))
747 d126363c Jim Regetz
  points(mean.circular(cx[1,], na.rm=TRUE), col="red")
748 8ceb6da6 Jim Regetz
  lines.circular(c(mean.circular(cx[1,], na.rm=TRUE), 0), c(0,-1),
749
      lwd=2, col="red")
750
  rose.diag(cx[2,], bins=1000, border="blue", ticks=FALSE, axes=FALSE)
751
  mtext(paste("(b)", round(yFromRow(sfd.bg, y2), 3), "degrees"))
752
  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),
753
      labels=c("N", "W", "S", "E"))
754 d126363c Jim Regetz
  points(mean.circular(cx[2,], na.rm=TRUE), col="red")
755 8ceb6da6 Jim Regetz
  lines.circular(c(mean.circular(cx[2,], na.rm=TRUE), 0), c(0,-1),
756
      lwd=2, col="red")
757
@
758
  \label{rose-diag-flowdir}
759
\end{figure}
760 0ffc32ec Jim Regetz
761
%
762 e6e6f64f Jim Regetz
% Figures: Correlations and RMSEs
763 0ffc32ec Jim Regetz
%
764 8ceb6da6 Jim Regetz
\begin{figure}
765 0ffc32ec Jim Regetz
  \centering
766 8ceb6da6 Jim Regetz
  \caption{Elevation associations between original and fused layers}
767
  \subfloat[Elevation correlations]{
768
    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../dem/elevation-assessment.pdf}
769
    \label{corr-elevation}
770 0ffc32ec Jim Regetz
  }\\
771 8ceb6da6 Jim Regetz
  \subfloat[Elevation RMSEs]{
772
    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../dem/elevation-assessment.pdf}
773
    \label{rmse-elevation}
774 0ffc32ec Jim Regetz
  }
775
\end{figure}
776 8ceb6da6 Jim Regetz
777
\begin{figure}
778 0ffc32ec Jim Regetz
  \centering
779 8ceb6da6 Jim Regetz
  \caption{Slope associations between original and fused layers}
780
  \subfloat[Slope correlations]{
781
    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../slope/slope-assessment.pdf}
782
    \label{corr-slope}
783 0ffc32ec Jim Regetz
  }\\
784 8ceb6da6 Jim Regetz
  \subfloat[Slope RMSEs]{
785
    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../slope/slope-assessment.pdf}
786
    \label{rmse-slope}
787 3c841968 Jim Regetz
  }
788
\end{figure}
789
790 8ceb6da6 Jim Regetz
\begin{figure}
791
  \centering
792
  \caption{Aspect associations between original and fused layers}
793
  \subfloat[Aspect circular correlations]{
794
    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../aspect/aspect-assessment.pdf}
795
    \label{corr-aspect}
796
  }\\
797
  \subfloat[Aspect RMSEs]{
798
    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../aspect/aspect-assessment.pdf}
799
    \label{rmse-aspect}
800
  }
801
\end{figure}
802 0ffc32ec Jim Regetz
803 8ceb6da6 Jim Regetz
\begin{figure}
804
  \centering
805
  \caption{Flow direction associations between original and fused layers}
806
  \subfloat[Flow direction correlations]{
807
    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../flow/flowdir-assessment.pdf}
808 e6e6f64f Jim Regetz
    \label{corr-flowdir}
809 8ceb6da6 Jim Regetz
  }\\
810
  \subfloat[Flow direction RMSEs]{
811
    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../flow/flowdir-assessment.pdf}
812 e6e6f64f Jim Regetz
    \label{rmse-flowdir}
813 8ceb6da6 Jim Regetz
  }
814
\end{figure}
815 3c841968 Jim Regetz
816
817 8711607d Jim Regetz
\clearpage
818
\appendix
819 e6e6f64f Jim Regetz
%-----------------------------------------------------------------------
820 8ceb6da6 Jim Regetz
\section{Code listings}
821 e6e6f64f Jim Regetz
%-----------------------------------------------------------------------
822 8ceb6da6 Jim Regetz
823 d126363c Jim Regetz
\lstinputlisting[language=bash, caption={GDAL commands for assembling
824
    and resampling SRTM and ASTER tiles into GeoTIFFs for use as inputs
825
    to the boundary correction routines described in this document.},
826
    label=code-gdalwarp]{../dem/boundary-assembly.sh}
827 8ceb6da6 Jim Regetz
\clearpage
828 d126363c Jim Regetz
829
\lstinputlisting[language=R, caption={R code implementing several
830
    SRTM-ASTER boundary corrections, including the Gaussian weighted
831
    average layer discussed in this document (see \texttt{OPTION 3}, as
832
    identified in code comments).},
833
    label=code-correct]{../dem/boundary-correction.R}
834
\clearpage
835
836
\lstinputlisting[language=bash, caption={GDAL commands for assembling
837
    boundary-corrected DEM components above and below 60\N, for each of
838
    several correction approaches implemented in Listing
839
    \ref{code-correct}. Note that multiresolution spline blending is
840
    treated separately (see Listing \ref{code-enblend}).},
841
    label=code-fuse]{../dem/boundary-fusion.sh}
842
\clearpage
843
844
\lstinputlisting[language=R, caption={R wrapper code for multiresolution
845
    spline of SRTM and ASTER.}, label=code-enblend]{../dem/enblend.R}
846
\clearpage
847
848 8ceb6da6 Jim Regetz
\lstinputlisting[language=bash, caption={GRASS code for computing
849 d126363c Jim Regetz
    flow directions using the various fused elevation rasters and their
850
    components DEMS as inputs.},
851
    label=code-flowdir]{../flow/flow-boundary.sh}
852 8ceb6da6 Jim Regetz
853 3c841968 Jim Regetz
\end{document}