% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-

%

% Using Gregor Gorjanc's bash script:

%   sweave -old=<pdfviewer> <filename.Rnw>

%

% Using my simple bash script wrapper:

%   sweave_jr <basename>

%

% Using the lower level commands, in two steps: 

%   R CMD Sweave <filename.Rnw>

%   pdflatex <filename.tex>

%

% To extract R code from within R:

%   Stangle("<filename.Rnw>", annotate=FALSE)

\usepackage{a4wide}

\usepackage{verbatim}

\setlist{noitemsep}

\lstset{

  commentstyle=\itshape\ttfamily,

  showstringspaces=false,

}

\newcommand{\N}{\ensuremath{^\circ\mbox{N}}\xspace}

\newcommand{\W}{\ensuremath{^\circ\mbox{W}}\xspace}

% add page numbers, but no headers

\thispagestyle{plain}

\title{SRTM/ASTER boundary analysis}

\date{Last update: 13 Jul 2011}

% don't number document sections

\setcounter{secnumdepth}{-1}. 

<<echo=FALSE,results=hide>>=

  # Gaussian weighted average parameters

  gwa.r <- 0.001

  gwa.50.pix <- round(sqrt(-log(0.5)/gwa.r))

  gwa.50.km <- gwa.50.pix * 90 / 1000

  gwa.01.pix <- round(sqrt(-log(0.01)/gwa.r))

  gwa.01.km <- gwa.01.pix * 90 / 1000

  # mean/median ASTER-SRTM elevation differences

  library(raster)

  demdir <- "/home/regetz/media/temp/terrain/dem/"

  d.srtm <- raster(file.path(demdir, "srtm_150below.tif"))

  d.aster <- raster(file.path(demdir, "aster_300straddle.tif"))

  d.delta.vals <- values(crop(d.aster, extent(d.srtm))) - values(d.srtm)

  delta.median <- median(d.delta.vals)

  delta.mean <- mean(d.delta.vals)

  delta.sd <- sd(d.delta.vals)

  delta.q <- quantile(d.delta.vals, c(0.1, 0.25, 0.75, 0.9))

@

\paragraph{Brief summary of findings}

   \item \ldots but with variability: standard deviation of

     $\mbox{ASTER}_i-\mbox{SRTM}_i$ is \Sexpr{round(delta.sd, 1)}

   \item ASTER also has numerous spurious spikes

     affecting slope/aspect?

  \end{itemize}

   \item (\emph{no further treatment in this document})

   \item leaves abrupt transition in derived slope and aspect

   \item unclear whether derived values of aspect and flow direction

     in the transition zone are acceptable

   \item also smooths transition slope and aspect

   \item unclear whether derived values of aspect and flow direction in

     the transition zone are acceptable

 \item Canada DEM itself has problems

  \begin{itemize}

   \item 60\N coincides with provincial boundaries; there are

     clear 60\N artifacts in this layer!

   \item other evident tiling artifacts too

  \end{itemize}

 \item Other comments

  \begin{itemize}

   \item N/S bias to aspect, flow direction computed on unprojected data

     at higher latitudes?

  \end{itemize}

\end{itemize}

\paragraph{To do (possibly)}

\begin{itemize}

 \item apply low pass filter to ASTER to reduce high frequency

 variation?

 \item apply algorithm to remove spikes (\ldots but maybe beyond scope?)

boundary between SRTM and ASTER, I focused on a narrow band along the

60\N boundary in Canada (Figure \ref{focal-area}). The latitudinal

extent of this band is 59.875\N - 60.125\N (i.e., 300 3" pixels

straddling 60\N), and the longitudinal extent is 136\W to 96\W (i.e.,

48,000 pixels wide). See Listing \ref{code-gdalwarp} for code. Within

this focal region, I generated latitudinal profiles of mean elevation,

slope, aspect, and flow direction, using the separate SRTM and ASTER

component DEMs themselves, using several different fused ASTER/SRTM DEMs

(see below), and using the Canadian Digital Elevation Data (CDED) as an

independent reference layer

(\url{http://www.geobase.ca/geobase/en/data/cded/index.html}). I also

then computed latitudinal correlations and RMSEs between the fused

layers and each of SRTM, ASTER, and CDED.

\paragraph{Elevation} This document includes latitudinal

characterizations of terrain values based on three different variants of

assessment of these here.

boundary artifacts (Figure \ref{blend-simple}).

(1983) method for blending overlapping images using multiresolution

4.0, \url{http://enblend.sourceforge.net}). Data preparation and

post-processing were handled in R (see Listing \ref{code-enblend}). As

presented here, the SRTM and ASTER inputs were prepared such that the

overlap zone is 75 latitudinal rows (6.75km) (Figure

\ref{blend-multires}).

weighted averaging such that the relative contribution of the SRTM

moving south away from 60\N (Equation \ref{eq-gaussian}).

\label{eq-gaussian}

&fused_{x,y} =

  \left\{

    \begin{array}{c l}

      ASTER_{x,y} & \mbox{above } 60^\circ \mbox{N} \\

      w_{x,y}ASTER_{x,y} + (1-w_{x,y})SRTM_{x,y} & \mbox{below } 60^\circ \mbox{N}

    \end{array}

  \right.\\

&\mbox{where }

  w_{x,y}=e^{-rD_{y}^{2}}

  \mbox{ and }

  D_{y} \mbox{ is the distance from } 60^\circ \mbox{N in units of pixels.} \nonumber

\end{eqnarray}

parameterized using $r$=\Sexpr{gwa.r}, producing equal weights for SRTM

and ASTER at a distance of $\sim$\Sexpr{round(gwa.50.km, 1)}km

(\Sexpr{gwa.50.pix} cells) south of the the boundary, and a relative

weight for ASTER of only 1\% by $\sim$\Sexpr{round(gwa.01.km, 1)}km

(\Sexpr{gwa.01.pix} cells) (Figure \ref{blend-gaussian}). See

\texttt{OPTION 3} in Listing \ref{code-correct}.

\subparagraph{Others not shown} I also experimented with some additional

fusion approaches, but have excluded them from further analysis in this

document.

\emph{Fused with exponential ramp north of 60\N.}

The first step was to take the pixel-wise difference between SRTM and

then added back into the ASTER values north of 60\N. This does a fine

job of eliminating the artificial shelf and thus the appearance of a

seam right along the 60\N boundary, but it does not address the abrupt

variability moving north across the boundary). Additionally, it

introduces vertical ``ridges'' running north from the boundary. These

arise because the calculated ramps are independent from one longitudinal

``column'' to the next, and thus any changes in the boundary difference

from one pixel to the next lead to adjacent ramps with different

inclines.

This involved first calculating the difference between SRTM and ASTER

everywhere south of 60\N, and then fitting a LOESS curve to these

differences using the actual ASTER elevation as a predictor. I then used

the fitted model to predict the ASTER-SRTM difference for each ASTER

cell north of 60\N, and added a declining fraction (based on a Gaussian

curve) of this difference to the corresponding ASTER elevation.

Conceptually, this amounts to applying an ASTER-predicted SRTM

correction to the ASTER elevation values, where the correction term has

a weight that declines to zero with increasing distance (north) away

from the boundary. However, this method alone didn't yield particularly

promising results in removing the 60\N seam itself, presumably because

adding a predicted correction at the boundary does not close the

SRTM-ASTER gap nearly as efficiently as do corrections based directly on

the observed SRTM vs ASTER elevation differences. I therefore haven't

useful in combination with one of the other methods.

For each of the three main fused DEM variants described above, slope was

    $ gdaldem slope -s 111120 <input_elevation> <output_slope>

\end{verbatim}

Note that the scale option used here is as recommended in the

\texttt{gdaldem} documentation:

\begin{quote}

  Lat/Long WGS84 projection), you can use scale=111120 if the vertical'

  units are meters}''

The output slope raster is in units of degrees.

\paragraph{Aspect}

As was the case with slope, aspect was calculated using

\texttt{gdaldem}:

\begin{verbatim}

    $ gdaldem aspect -s 111120 <input_elevation> <output_aspect>

\end{verbatim}

The output aspect raster values indicate angular direction in units of

degrees, with 0=North and proceeding clockwise.

\paragraph{Flow direction}

Flow direction was calculated using the GRASS (GRASS GIS 6.4.1)

\texttt{r.terraflow} module; see code listing \ref{code-flowdir}.

Because of a $\sim$30k ($2^{15}$) limit to the input raster dimension

size in the pre-built GRASS \texttt{r.terraflow} module I used, this

analysis was restricted to a smaller longitudinal subset of the data,

spanning 125\W to 100\W.

The default flow direction output of this module is encoded so as to

indicate \emph{all} downslope neighbors, also known as the Multiple Flow

Direction (MFD) model. However, to simplify post-processing and

summarization, the results here are based on an alternative Single Flow

Direction (SFD, \emph{a.k.a.} D8) model, which indicates the neighbor

associated with the steepest downslope gradient. Note that this is

equivalent to what ArcGIS GRID \texttt{flowaccumulation} command does. I

then recoded the output raster to use the same azimuth directions used

by \texttt{gdaldem aspect}, as described above for aspect.

\section{Latitudinal mean terrain profiles}

%-----------------------------------------------------------------------

(Figure \ref{mean-elevation}). SRTM tends to be highest, ASTER is

between SRTM and ASTER is fairly consistent not only across latitudes,

but also across elevations (Figure \ref{aster-srtm}). The overall median

difference between ASTER and SRTM (i.e., considering

$\mbox{ASTER}_i-\mbox{SRTM}_i$ for all pixels $i$ where the two DEMs

co-occur) is \Sexpr{delta.median} meters, with a mean of

a few meters) across the observed range of elevations (Figure

consistent across latitudes and across elevation zones, additional

variation is evident at the pixel level. Again focusing on the

pixel-wise differences, they appear to be symmetrically distributed

about the mean with a standard deviation of \Sexpr{round(delta.sd, 1)}

and quartiles ranging from \Sexpr{delta.q["25%"]} to

\Sexpr{delta.q["75%"]} meters; ASTER elevations are actually greater

than SRTM for \Sexpr{round(100*mean(d.delta.vals>0))}\% of pixels (see

Figure \ref{aster-srtm-scatter}). Thus, although adding a constant

offset of \Sexpr{-delta.median} meters to the ASTER DEM would clearly

\ref{aster-srtm-scatter} also highlights the existence of several

obviously spurious ASTER spikes of >1000m; although not shown here,

these tend to occur in small clumps of pixels, perhaps corresponding to

false elevation readings associated with clouds?

Not surprisingly, simple fusion produces an artificial $\sim$12m cliff

in the mean elevation profile (Figure \ref{mean-elevation}). At least in

terms of mean elevation, this artifact is completely removed by both the

although ultimately this would depend on the chosen zone of overlap and

on the exact parameterization of the weighting function.

\paragraph{Slope} The mean ASTER slope is uniformly steeper than the

mean SRTM slope at all latitudes in the area of overlap, by nearly 1

degree (Figure \ref{mean-slope}). However, the shape of the profile

itself is nearly identical between the two. Although this may partly

reflect inherent SRTM vs ASTER differences, my guess is that CGIAR

\textbf{\color{red}[todo: check!]}

ASTER (presumably because it is at least partially derived from

contour-based data \textbf{\color{red}[todo: check!]}). Moreover, this

figure makes it clear that CDED has some major artifacts at

provincial boundaries across the entirety of western Canada) means we

The simple fused layer exhibits a dramatic spike in slope at the

is smoothed out in the latter. This likely reflects the fact that the

multiresolution spline effectively uses a very narrow transition zone

for stitching together high frequency components of the input images,

for the shift in mean slope. 

\paragraph{Aspect}

For the purposes of latitudinal profiles, aspect values were summarized

using a circular mean (Equation \ref{eq-cmean}).

\begin{equation}

 \bar{x} = \mathrm{atan2}

   \left(

    \sum_{i=1}^{n}\frac{\sin(x_i)}{n},

    \sum_{i=1}^{n}\frac{\cos(x_i)}{n}

   \right )

\label{eq-cmean}

\end{equation}

where $x_i$ is the aspect value (in radians) of pixel $i$.

As indicated in Figure \ref{mean-aspect}, the circular mean aspect

the area of overlap, and mean aspect values calculated on CDED

latitudes are directed either nearly north or nearly south, though

almost always with a slight eastward rather than westward inclination.

In fact, there appears to be a general bias towards aspect values

orienting along the north-south axis, as is especially apparent in the

rose diagrams of Figure \ref{rose-diag-aspect}. I suspect this is an

artifact of our use of unprojected data, especially at these high

latitudes. Because the unprojected raster is effectively 'stretched'

east and west relative to the actual topography, elevational gradients

along the east-west axis are artificially flattened, and the direction

of dominant gradient is more likely to be along the north-south axis.

Upon further reflection, it's not clear whether these patterns of mean

aspect are particular useful diagnostics, as they seems to be sensitive

to subtle variations in the data. Referring again to Figure

\ref{rose-diag-aspect}, note how the mean direction flips from nearly

north to nearly south between the two latitudes, even though the

distributions of pixel-wise aspect values are nearly indistinguishable

by eye.

In any case, not surprisingly, the simple fused layer matches the SRTM

aspect values south of 60\N and the ASTER aspect values north of 60\N; at

the immediate boundary, the mean aspect is northward, as one would

expect in the presence of a cliff artifact at the seam.

Interestingly, the aspect layers derived from the two blended DEMs

respective fusion zones. This pattern is visually obvious at latitudes

between 59.95\N and 60\N in the bottom two panels of Figure

\ref{mean-aspect}. This is almost certainly a signal of the blending of

the lower elevation ASTER to the north with higher elevation SRTM to the

throughout this zone. 

northward at all latitudes, and SRTM-derived flow is nearly northward at

all but a few latitudes (Figure \ref{mean-flowdir}). This seems

reasonable considering that most pixels in this Canada test region fall

in the Arctic drainage. For unknown reasons, CDED produces southward

mean flow direction at numerous latitudes, and generally seems to have a

slightly more eastward tendency. As was the case with aspect, note that

a general north-south bias is evident (Figure \ref{rose-diag-flowdir}),

again likely due to use of an unprojected raster at high latitudes.

The various fused layer profiles look as one would expect (Figure

\ref{mean-flowdir}), although the overall lack of latitudinal

variability in mean flow direction in SRTM, ASTER, and all three derived

layers makes it hard to say much more than that. 

%-----------------------------------------------------------------------

\section{Informal correlation analysis}

%-----------------------------------------------------------------------

\paragraph{Elevation}

>0.999, and RMSEs are on the order of 10-15 meters (Figures

\ref{corr-elevation} and \ref{rmse-elevation}). Spikes (downward for

correlation, upward for RMSE) occur at some latitudes, quite possibly

associated with the observed extreme spikes in the ASTER DEM itself. As

expected, the multiresolution spline and Gaussian weighted average both

produce layers that gradually become less similar to ASTER and more

similar to SRTM moving south from 60\N, but in slightly different ways.

This gradual transition is less evident when considering associations

with CDED (bottom panels of Figures \ref{corr-elevation} and

\paragraph{Slope}

The patterns of slope layer similarity are much like those described

above for elevation, although the correlations are somewhat lower

($\sim$0.94 between SRTM and ASTER (Figure \ref{corr-slope})). RMSEs

between SRTM and ASTER are approximately 2 at all latitudes where the

data co-occur (Figure \ref{rmse-slope}). Another difference, echoing a

pattern previously noted in the profiles of mean slope itself, is that

Gaussian weighted averaging produces a layer that exhibits a gradual

transition from ASTER to SRTM, whereas the multiresolution spline yields

an abrupt transition. Not surprisingly, the simple fused layer is even

and \ref{rmse-slope}.

\paragraph{Aspect}

Because aspect values are on a circular scale, I calculated modified

versions of the Pearson correlation coefficient using the

\textbf{circular} R package function \texttt{cor.circular}. I believe

this implements the formula described by Jammalamadaka \& Sarma (1988):

\begin{equation}

r = \frac

  {\sum_{i=1}^{n} \sin(x_i-\bar{x}) \sin(y_i-\bar{y})}

  {\sqrt{\sum_{i=1}^{n} \sin(x_i-\bar{x})^2}

   \sqrt{\sum_{i=1}^{n} \sin(y_i-\bar{y})^2}}

\label{eq-ccor}

\end{equation}

where $x_i$ and $y_i$ are aspect values (in radians) for pixel $i$, and

$\bar{x}$ and $\bar{y}$ are circular means calculated as in Equation

\ref{eq-cmean}.

properties analogous to computation of the circular correlation, but

instead just imposed a simple correction whereby all pairwise

differences were computed using the shorter of the two paths around the

compass wheel (Equation \ref{eq-circ-rmse}). For example, the difference

between 0$^\circ$ and 150$^\circ$ is 150$^\circ$, but the difference

between 0$^\circ$ and 250$^\circ$ is 110$^\circ$.

\begin{eqnarray}

\label{eq-circ-rmse}

&RMSE = \sqrt{\frac{\sum_{i=1}^{n} (\Delta_i^2)}{n}}\\ \nonumber

&\mbox{where }

  \Delta_i = \mbox{argmin} \left (|x_i-y_i|, 360-|x_i - y_i| \right )

\end{eqnarray}

Circular correlations between SRTM and ASTER were surprisingly low,

typically hovering around 0.5, but dipping down towards zero at numerous

latitudes (Figure \ref{corr-flowdir}). The corresponding RMSE is close

to 70 at all latitudes, surprisingly high considering that the maximum

difference between any two pixels is 180$^\circ$ (Figure

For both of the blended layers (multiresolution spline and Gaussian

weighted average), aspect values calculated in the zone of ASTER/SRTM

overlap appear to be \emph{less} similar to the component DEMs than

those to data sources are to each other. As evident in the upper middle

and upper right panels of Figure \ref{corr-flowdir}, correlations

between fused aspect and aspect based on the original SRTM and ASTER

images are substantially \emph{negative} for many latitudes in the zone

of overlap. I don't have a good feel for what's going on here, although

that I don't have a good feel for. Note that the latitudinal profile of

aspect RMSEs is less worrisome (upper middle and upper right panels of

Figure \ref{rmse-flowdir}) and more closely resembles the profile of

slope RMSEs discussed above, particularly as regards the more gradual

multiresolution spline.

\paragraph{Flow direction}

As with aspect, circular correlation coefficients and adjusted RMSEs

around 0.40 but spiking negative at a handful of latitudes (Figure

correlation profiles are fairly well-behaved for both blended layers,

and don't show the same odd behavior as was the case for aspect. Again

direction profiles echo this pattern (Figure \ref{rmse-flowdir}), and

indeed look almost the same as those computed using aspect (Figure

\ref{rmse-aspect}).

%-----------------------------------------------------------------------

% FIGURES

%-----------------------------------------------------------------------

\clearpage

%

% Figure: Boundary analysis region

%

\begin{figure}[h]

  that for flow direction, analysis was restricted to a smaller

  longitudinal span (125\W to 100\W).}

<<echo=FALSE,results=hide,fig=TRUE,height=5,width=7>>=

  par(omi=c(0,0,0,0))

  library(raster)

  library(maps)

  demdir <- "/home/regetz/media/temp/terrain/dem/"

  dem <- raster(file.path(demdir, "fused_300straddle.tif"))

  map("world", xlim=c(-140,-70), ylim=c(40, 70), fill=TRUE, col="grey",

      mar=c(4,0,0,0))

  grid(col="darkgray")

  axis(1)

  axis(2)

  box()

  plot(extent(dem), col="red", add=TRUE)

  arrows(-110, 57, -113, 59.7, col="red", length=0.1)

  text(-110, 57, labels="focal region", col="red", font=3, pos=1)

@

  \label{focal-area}

\end{figure}

%

% Figure: Fusion methods

%

\begin{figure}

  \caption{SRTM/ASTER DEM fusion configurations}

  \centering

  \subfloat[Simple fusion]{

    \label{blend-simple}

<<echo=FALSE,results=hide,fig=TRUE, height=3, width=6>>=

  par(omi=c(0,0,0,0), mar=c(4,4,0,2)+0.1)

  plot(0, xlim=c(-150, 150), ylim=c(0, 1), type="n", xaxt="n", yaxt="n",

      bty="n", xlab="Latitude", ylab=NA)

  rect(0, 0, 150, 0.75, col="red", density=5, angle=45)

  rect(-150, 0, 0, 0.75, col="blue", density=5, angle=135)

#  axis(1, at=seq(-120, 120, by=60), labels=seq(59.9, 60.1, by=0.05))

  axis(1, at=seq(-150, 150, by=75), labels=c(59.875, NA, 60, NA, 60.125))

  text(-150, 0.9, labels="SRTM", pos=4, col="blue")

  text(150, 0.9, labels="ASTER", pos=2, col="red")

@

  }\\

  \subfloat[Multiresolution spline]{

    \label{blend-multires}

<<echo=FALSE,results=hide,fig=TRUE, height=3, width=6>>=

  par(omi=c(0,0,0,0), mar=c(4,4,0,2)+0.1)

  plot(0, xlim=c(-150, 150), ylim=c(0, 1), type="n", xaxt="n", yaxt="n",

      bty="n", xlab="Latitude", ylab=NA)

  rect(-75, 0, 0, 0.75, col="lightgray")

  rect(-75, 0, 150, 0.75, col="red", density=5, angle=45)

  rect(-150, 0, 0, 0.75, col="blue", density=5, angle=135)

  #axis(1, at=seq(-150, 150, by=60), labels=seq(59.875, 60.125, by=0.05))

  #axis(1, at=0, labels=60, col="red", cex=0.85)

  #axis(1, at=seq(-120, 120, by=60), labels=seq(59.9, 60.1, by=0.05))

  axis(1, at=seq(-150, 150, by=75), labels=c(59.875, NA, 60, NA, 60.125))

  text(-150, 0.9, labels="SRTM", pos=4, col="blue")

  text(-37.5, 0.8, labels="overlap", pos=3, font=3)

  text(150, 0.9, labels="ASTER", pos=2, col="red")

@

  }\\

  \subfloat[Gaussian weighted average]{

    \label{blend-gaussian}

<<echo=FALSE,results=hide,fig=TRUE, height=3, width=6>>=

  par(omi=c(0,0,0,0), mar=c(4,4,1,2)+0.1)

  curve(exp(-0.001*x^2), from=-150, to=0, xlim=c(-150, 150), col="red",

      xaxt="n", bty="n", xlab="Latitude", ylab="Weighting")

  curve(1+0*x, from=0, to=150, add=TRUE, col="red")

  curve(1-exp(-0.001*x^2), from=-150, to=0, add=TRUE, col="blue")

  #axis(1, at=seq(-150, 150, by=60), labels=seq(59.875, 60.125, by=0.05))

  #axis(1, at=0, labels=60, col="red", cex=0.85)

  #axis(1, at=seq(-120, 120, by=60), labels=seq(59.9, 60.1, by=0.05))

  axis(1, at=seq(-150, 150, by=75), labels=c(59.875, NA, 60, NA, 60.125))

  text(-100, 0.6, labels="gaussian\nblend")

  text(-140, 0.9, labels="SRTM", col="blue", pos=4)

  text(-140, 0.1, labels="ASTER", col="red", pos=4)

@

  }

\end{figure}

%

% Figure: Mean latitudinal profiles

%

\begin{figure}

  \centering

  \caption{Mean elevation (m)}

    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../dem/elevation-assessment.pdf}

    \label{mean-elevation}

\end{figure}

\begin{figure}

  \centering

  \caption{Mean slope (degrees)}

    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../slope/slope-assessment.pdf}

    \label{mean-slope}

\end{figure}

\begin{figure}

  \centering

  \caption{Circular mean aspect (azimuth)}

    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../aspect/aspect-assessment.pdf}

    \label{mean-aspect}

\end{figure}

\begin{figure}

  \centering

  \caption{Circular mean flow direction}

    \includegraphics[page=1, trim=1in 1in 1in 1in, width=\linewidth]{../flow/flowdir-assessment.pdf}

    \label{mean-flowdir}

\end{figure}

%

% Figure: ASTER vs SRTM elevation patterns

%

\begin{figure}

  \caption{ASTER vs SRTM elevation comparisons}

  \centering

  \subfloat[Boxplots of the arithmetic difference in elevations (ASTER -

  numbers above each whisker indicates how many thousands of pixels are

  included in the corresponding summary. Dashed red line indicates the

  median difference across all pixels south of 60\N.]{

    \includegraphics[width=\linewidth, trim=0 0 0 0.5in, clip=TRUE]{../dem/aster-srtm-bins.png}

  }\\

  the 1:1 diagonal, and the parallel red line is offset lower by the

  }

  \label{aster-srtm}

\end{figure}

%

% Figures: Aspect rose diagrams

%

  \caption{Binned frequency distributions of aspect within two selected

    latitudinal strips: (a) within the SRTM-ASTER blend zone, and (b)

    south of the overlap zone. Aspect values here are based on the DEM

    produced via Gaussian weighted averaging, but similar patterns are

    evident using the other DEMs. Red spoke lines indicate the circular

    mean direction across all pixels at the associated latitude.}

  \centering

<<echo=FALSE,results=hide,fig=TRUE,height=5,width=7>>=

  library(circular)

  aspdir <- "/home/regetz/media/temp/terrain/aspect"

  a.bg <- raster(file.path(aspdir, "fused_300straddle_blendgau_a.tif"))

  par(mfrow=c(1,2), mar=c(0,0,4,0))

  y1 <- 200

  y2 <- 220

  cx <- circular(as.matrix(a.bg)[c(y1,y2),], units="degrees",

      rotation="clock", zero=pi/2)

  rose.diag(cx[1,], bins=8, axes=FALSE)

  mtext(paste("(a)", round(yFromRow(a.bg, y1), 3), "degrees"))

  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),

      labels=c("N", "W", "S", "E"))

  points(mean.circular(cx[1,], na.rm=TRUE), col="red")

  lines.circular(c(mean.circular(cx[1,], na.rm=TRUE), 0), c(0,-1),

      lwd=2, col="red")

  rose.diag(cx[2,], bins=8, axes=FALSE)

  mtext(paste("(b)", round(yFromRow(a.bg, y2), 3), "degrees"))

  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),

      labels=c("N", "W", "S", "E"))

  points(mean.circular(cx[2,], na.rm=TRUE), col="red")

  lines.circular(c(mean.circular(cx[2,], na.rm=TRUE), 0), c(0,-1),

      lwd=2, col="red")

@

  \label{rose-diag-aspect}

\end{figure}

% Figures: Flow direction rose diagrams

%

  \caption{Relative frequencies of D8 flow directions (indicated by

    relative heights of the blue spokes) within two selected

    latitudinal bands: (a) within the SRTM-ASTER blend zone, and (b)

    south of the overlap zone. Flow directions here are based on the DEM

    produced via Gaussian weighted averaging. Red spoke lines indicate

    the circular mean flow direction across all pixels at the associated

    latitude.}

  \centering

<<echo=FALSE,results=hide,fig=TRUE,height=5,width=7>>=

  library(circular)

  flowdir <- "/home/regetz/media/temp/terrain/flow"

  # create function to recode terraflow SFD values into degrees, with

  # 0=North and proceeding clockwise (this matches gdaldem's default

  # azimuth output for aspect calculation)

  recode <- function(r) {

      v <- values(r)

      v[v==0] <- NA

      v[v==1] <- 90  ## east

      v[v==2] <- 135

      v[v==4] <- 180  ## south

      v[v==8] <- 225

      v[v==16] <- 270  ## west

      v[v==32] <- 315

      v[v==64] <- 0  ## north

      v[v==128] <- 45

      r[] <- v

      return(r)

  }

  sfd.bg <- recode(raster(file.path(flowdir,

      "fused_300straddle_blendgau_sfd.tif")))

  par(mfrow=c(1,2), mar=c(0,0,4,0))

  y1 <- 200

  y2 <- 220

  cx <- circular(as.matrix(sfd.bg)[c(y1,y2),], units="degrees",

      rotation="clock", zero=pi/2)

  rose.diag(cx[1,], bins=1000, border="blue", ticks=FALSE, axes=FALSE)

  mtext(paste("(a)", round(yFromRow(sfd.bg, y1), 3), "degrees"))

  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),

      labels=c("N", "W", "S", "E"))

      lwd=2, col="red")

  rose.diag(cx[2,], bins=1000, border="blue", ticks=FALSE, axes=FALSE)

  mtext(paste("(b)", round(yFromRow(sfd.bg, y2), 3), "degrees"))

  axis.circular(at=circular(c(pi/2, pi, 3*pi/2, 2*pi)),

      labels=c("N", "W", "S", "E"))

      lwd=2, col="red")

@

  \label{rose-diag-flowdir}

\end{figure}

%

  \subfloat[Elevation correlations]{

    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../dem/elevation-assessment.pdf}

    \label{corr-elevation}

    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../dem/elevation-assessment.pdf}

    \label{rmse-elevation}

\end{figure}

\begin{figure}

  \subfloat[Slope correlations]{

    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../slope/slope-assessment.pdf}

    \label{corr-slope}

    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../slope/slope-assessment.pdf}

    \label{rmse-slope}

\end{figure}

  \centering

  \caption{Aspect associations between original and fused layers}

  \subfloat[Aspect circular correlations]{

    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../aspect/aspect-assessment.pdf}

    \label{corr-aspect}

  }\\

  \subfloat[Aspect RMSEs]{

    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../aspect/aspect-assessment.pdf}

    \label{rmse-aspect}

  }

\end{figure}

  \centering

  \caption{Flow direction associations between original and fused layers}

  \subfloat[Flow direction correlations]{

    \includegraphics[page=3, trim=1in 1in 1in 1in, width=\linewidth]{../flow/flowdir-assessment.pdf}

  \subfloat[Flow direction RMSEs]{

    \includegraphics[page=2, trim=1in 1in 1in 1in, width=\linewidth]{../flow/flowdir-assessment.pdf}

\end{figure}

\appendix

    and resampling SRTM and ASTER tiles into GeoTIFFs for use as inputs

    to the boundary correction routines described in this document.},

    label=code-gdalwarp]{../dem/boundary-assembly.sh}

\lstinputlisting[language=R, caption={R code implementing several

    SRTM-ASTER boundary corrections, including the Gaussian weighted

    average layer discussed in this document (see \texttt{OPTION 3}, as

    identified in code comments).},

    label=code-correct]{../dem/boundary-correction.R}

\clearpage

\lstinputlisting[language=bash, caption={GDAL commands for assembling

    boundary-corrected DEM components above and below 60\N, for each of

    several correction approaches implemented in Listing

    \ref{code-correct}. Note that multiresolution spline blending is

    treated separately (see Listing \ref{code-enblend}).},

    label=code-fuse]{../dem/boundary-fusion.sh}

\clearpage

\lstinputlisting[language=R, caption={R wrapper code for multiresolution

    spline of SRTM and ASTER.}, label=code-enblend]{../dem/enblend.R}

\clearpage

    components DEMS as inputs.},

    label=code-flowdir]{../flow/flow-boundary.sh}