% -*- 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)

\documentclass{article}

\usepackage[font=small,labelfont=bf]{caption}

\usepackage[font=normalsize,singlelinecheck=false]{subfig}

\usepackage{enumitem}

\usepackage[utf8]{inputenc}

\usepackage{a4wide}

\usepackage{verbatim}

\usepackage{listings}

\usepackage[colorlinks=true,urlcolor=red]{hyperref}

\usepackage{color}

\usepackage{xspace}

% shrink list item spacing

\setlist{noitemsep}

% set some listings options

\lstset{

  basicstyle=\footnotesize\ttfamily,

  stringstyle=\ttfamily,

  commentstyle=\itshape\ttfamily,

  showstringspaces=false,

}

% define degree command for convenience

\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}

\author{Jim Regetz, NCEAS}

\date{Last update: 13 Jul 2011}

\begin{document}

\maketitle

% don't number document sections

\setcounter{secnumdepth}{-1}. 

% first compute some values to use in the text

<<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}

\begin{itemize}

 \item SRTM vs ASTER differences

  \begin{itemize}

   \item ASTER is systematically lower, by 12 meters in the median case

   \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

   \item ASTER has more high-frequency variability (``texture''),

     affecting slope/aspect?

  \end{itemize}

 \item Fusion via northward exponential rampdown of boundary delta

  \begin{itemize}

   \item eliminates elevation cliff at 60\N

   \item leaves abrupt transition in SRTM/ASTER textural differences

   \item introduces north-south ridging artifacts

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

  \end{itemize}

 \item Fusion via multiresolution spline

  \begin{itemize}

   \item eliminates elevation cliff at 60\N

   \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

  \end{itemize}

 \item Fusion via Gaussian weighted average of SRTM/ASTER

  \begin{itemize}

   \item eliminates elevation cliff at 60\N

   \item also smooths transition slope and aspect

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

     the transition zone are acceptable

  \end{itemize}

 \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 add constant offset of \Sexpr{-delta.median}m to ASTER

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

 variation?

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

\end{itemize}

\clearpage

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

\section{Terrain layer production methodology}

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

For the purposes of assessing artifacts associated with the northern

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

a fused 3" ASTER/SRTM DEM. I also explored (and briefly describe below) two

additional approaches to fusing the layers, but do not include further

assessment of these here.

\subparagraph{Simple fusion} Naive concatenation of SRTM below 60\N with

ASTER above 60\N, without applying any modifications to deal with

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

\subparagraph{Multiresolution spline} Application of Burt \& Adelson's

(1983) method for blending overlapping images using multiresolution

splines, as implemented in the \emph{Enblend} software package (version

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}).

\subparagraph{Gaussian weighted average} Blend of the two layers using

weighted averaging such that the relative contribution of the SRTM

elevation is zero at 60\N, and increases as a function of distance

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

\begin{eqnarray}

\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}

For the assessment presented here, the weighting function function was

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

ASTER in the row immediately below 60\N (i.e., the northernmost extent of

SRTM). An exponentially declining fraction of this difference was then

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

transition in texture (i.e., the sudden appearance of high frequency

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.

\emph{Simple LOESS predictive model.}

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

pursued this any further, although it (or something like it) may prove

useful in combination with one of the other methods.

\paragraph{Slope}

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

calculated using \texttt{gdaldem} (GDAL 1.8.0, released 2011/01/12):

\begin{verbatim}

    $ 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}

  ``\emph{If the horizontal unit of the source DEM is degrees (e.g

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

  units are meters}''

\end{quote}

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}

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

\paragraph{Elevation} SRTM, ASTER, and CDED all share a very

similarly shaped mean elevation profile, but with differing heights

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

lowest, and CDED is intermediate. The magnitude of average difference

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

\Sexpr{round(delta.mean, 2)} meters, and this more or less holds (within

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

\ref{aster-srtm-bins}). However, while this average offset is broadly

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

center it with respect to the SRTM (at least in the Canada focal

region), appreciable differences would remain. Figure

\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

multiresolution spline and Gaussian weighted average methods. The

transition is, to the eye, slightly smoother in the former case,

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

post-processing of the particular SRTM product we're using has removed

some of the high frequency ``noise'' that remains in ASTER?

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

Note that the CDED tends to be flatter than both SRTM and

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

regular intervals. The spike especially at 60\N (which coincides with

provincial boundaries across the entirety of western Canada) means we

probably need to scuttle our plans to use this DEM as a formal reference

dataset for boundary analysis.

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

immediate 60\N boundary, undoubtedly associated with the artificial

elevation cliff. This artifact is eliminated by both the multiresolution

spline and Gaussian weighting. However, the former exhibits a sudden step

change in slope in the SRTM-ASTER overlap region, whereas the transition

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,

and it seems likely that these are precisely the features responsible

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

values of SRTM and ASTER are generally similar across all latitudes in

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

are similar at most but not all latitudes. The mean values at nearly all

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

(multiresolution spline and Gaussian weighted average) exhibit a

consistent mean northward inclination at all latitudes in their

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

south, introducing a north-facing tilt (however slight) to the data

throughout this zone. 

\paragraph{Flow direction} With the exception of edge effects at the

margins of the input rasters, mean ASTER-derived flow is nearly

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}

Pearson correlations between SRTM and ASTER are quite high, typically

>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

\ref{rmse-elevation}), which in general is less correlated with SRTM,

and even less with ASTER, than those two layers are with each other.

\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

worse, producing not only a sudden transition but also aberrant values

at the fusion seam itself; note downward (upward) spikes in correlation

(RMSE) at 60\N in the first column of plots in Figures \ref{corr-slope}

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}.

To calculate RMSEs for aspect, I did not attempt to use trigonometric

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

\ref{rmse-flowdir}). Comparison of SRTM and ASTER with CDED

yields similar patterns.

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

it may just involve properties of the circular correlation statistic

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

transition in case of Gaussian weighted averaging compared to the

multiresolution spline.

\paragraph{Flow direction}

As with aspect, circular correlation coefficients and adjusted RMSEs

were calculated for each latitudinal band. Flow direction correlations

between SRTM and ASTER are slightly lower than for aspect, typically

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

\ref{corr-flowdir}). RMSEs hover consistently around $\sim$75, slightly

lower than was the case with aspect (Figure \ref{rmse-flowdir}).

Aside from an expected flow artifact at the southern image edge, and an

unexpected (and as-yet unexplained) negative spike at $\sim$59.95\N, the

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

it is clear that the multiresolution spline results in a much more abrupt

transition than does the Gaussian weighted average. The RMSE flow

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]

  \caption{Focal area used as the basis for boundary assessment. Note

  that for flow direction, analysis was restricted to a smaller

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

  \centering

<<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 -

  SRTM), summarized across a series of SRTM elevation bins. Boxes

  include the median (horizontal band), 1st and 3rd quartiles (box

  extents), and $\pm$1.5$\times$IQR (whiskers); outliers not shown. Gray

  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}

    \label{aster-srtm-bins}

  }\\

  \subfloat[Pixel-wise plot of SRTM vs ASTER for all values in the 150

  latitudinal rows of overlap south of 60\N. Dashed blue line indicates

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

  observed median difference between the two DEMs (\Sexpr{-delta.median}m).

  Inset histogram shows distribution of differences, excluding absolute

  differences >60m.]{

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

    \label{aster-srtm-scatter}

  }

  \label{aster-srtm}

\end{figure}

%

% Figures: Aspect rose diagrams

%

\begin{figure}[h!]

  \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

%

\begin{figure}[h!]

  \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"))

  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=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"))

  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-flowdir}

\end{figure}

%

% Figures: Correlations and RMSEs

%

\begin{figure}

  \centering

  \caption{Elevation associations between original and fused layers}

  \subfloat[Elevation correlations]{

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

    \label{corr-elevation}

  }\\

  \subfloat[Elevation RMSEs]{

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

    \label{rmse-elevation}

  }

\end{figure}

\begin{figure}

  \centering

  \caption{Slope associations between original and fused layers}

  \subfloat[Slope correlations]{

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

    \label{corr-slope}

  }\\

  \subfloat[Slope RMSEs]{

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

    \label{rmse-slope}

  }

\end{figure}

\begin{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}

\begin{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}

    \label{corr-flowdir}

  }\\

  \subfloat[Flow direction RMSEs]{

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

    \label{rmse-flowdir}

  }

\end{figure}

\clearpage

\appendix

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

\section{Code listings}

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

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

    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}

\clearpage

\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

\lstinputlisting[language=bash, caption={GRASS code for computing

    flow directions using the various fused elevation rasters and their

    components DEMS as inputs.},

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

\end{document}