# R code to plot latitudinal profiles of (circular) mean aspect, along

# with both RMSE and (circular) correlation coefficients comparing fused

# layers with both the raw ASTER and with the Canada DEM

#

# Aspect layers were generated from the respect DEMs in this way:

#   $ gdaldem aspect -s 111120 <layer>.tif <layer>_a.tif

# ...where the default azimuthal behavior produces output values ranging

# from 0-360 where 0 is north, and proceeding clockwise

#

# For exploratory plotting, note the following (uses 'circular'

# package):

#  > cx <- circular(as.matrix(a.bg)[151,], units="degrees",

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

#  > rose.diag(cx, bins=8)

#  > points(mean.circular(cx, na.rm=TRUE), col="red")

#

# Jim Regetz

# NCEAS

library(raster)

library(circular)

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

# load aspect rasters

a.aster <- raster(file.path(datadir, "aster_300straddle_a.tif"))

a.srtm <- raster(file.path(datadir, "srtm_150below_a.tif"))

a.uncor <- raster(file.path(datadir, "fused_300straddle_a.tif"))

a.enblend <- raster(file.path(datadir, "fused_300straddle_enblend_a.tif"))

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

a.can <- raster(file.path(datadir, "cdem_300straddle_a.tif"))

# extract raster latitudes for later

lats300 <- yFromRow(a.aster, 1:nrow(a.aster))

lats150 <- yFromRow(a.srtm, 1:nrow(a.srtm))

# initialize output pdf device driver

pdf("aspect-assessment.pdf", height=8, width=11.5)

#

# plot latitudinal profiles of mean aspect

#

# simple helper function to calculate row-wise means using circular

# mean, patterned after circ.mean in the CircStats package

rowMeansC <- function(r1, na.rm=TRUE) {

    m1 <- as.matrix(r1)

    m1[] <- (m1 * pi)/180

    sinr <- rowSums(sin(m1), na.rm=na.rm)

    cosr <- rowSums(cos(m1), na.rm=na.rm)

    cmeans <- atan2(sinr, cosr)

    (cmeans * 180)/pi

}

par(mfrow=c(2,2), omi=c(1,1,1,1))

ylim <- c(-180, 180)

plot(lats300, rowMeansC(a.can), type="l", yaxt="n",

    xlab="Latitude", ylab="Mean aspect", ylim=ylim)

axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Original DEMs")

lines(lats300, rowMeansC(a.aster), col="blue")

lines(lats150, rowMeansC(a.srtm), col="red")

legend("bottomright", legend=c("ASTER", "SRTM", "CDED"), col=c("blue",

    "red", "black"), lty=c(1, 1), bty="n")

abline(v=60, col="red", lty=2)

mtext(expression(paste("Latitudinal profiles of mean aspect (",

    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)

plot(lats300, rowMeansC(a.uncor), type="l", yaxt="n",

    xlab="Latitude", ylab="Mean aspect", ylim=ylim)

axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="simple fused")

abline(v=60, col="red", lty=2)

plot(lats300, rowMeansC(a.enblend), type="l", yaxt="n",

    xlab="Latitude", ylab="Mean aspect", ylim=ylim)

axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="multires spline")

abline(v=60, col="red", lty=2)

plot(lats300, rowMeansC(a.bg), type="l", yaxt="n",

    xlab="Latitude", ylab="Mean aspect", ylim=ylim)

axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend")

abline(v=60, col="red", lty=2)

#

# plot latitudinal profiles of RMSE

#

# simple helper function to calculate row-wise RMSEs, accounting for the

# fact that aspect values are circular (0-360), so the difference

# between e.g. 5 and 355 should only be 10

rmse <- function(r1, r2, na.rm=TRUE, use) {

    diffs <- abs(as.matrix(r1) - as.matrix(r2))

    if (!missing(use)) diffs[!use] <- NA

    diffs[] <- ifelse(diffs>180, 360-diffs, diffs)

    sqrt(rowMeans(diffs^2, na.rm=na.rm))

}

par(mfrow=c(2,3), omi=c(1,1,1,1))

ylim <- c(0, 100)

# ...with respect to ASTER

plot(lats300, rmse(a.uncor, a.aster), type="l", xlab="Latitude",

    ylab="RMSE", ylim=ylim)

lines(lats150, rmse(crop(a.uncor, extent(a.srtm)), a.srtm), col="blue")

legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),

    lty=c(1, 1), bty="n")

text(min(lats300), max(ylim)-5, pos=4, font=3, labels="simple fused")

abline(v=60, col="red", lty=2)

mtext(expression(paste(

    "Aspect discrepancies with respect to separate ASTER/SRTM components (",

    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)

plot(lats300, rmse(a.enblend, a.aster), type="l", xlab="Latitude",

    ylab="RMSE", ylim=ylim)

lines(lats150, rmse(crop(a.enblend, extent(a.srtm)), a.srtm), col="blue")

legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),

    lty=c(1, 1), bty="n")

text(min(lats300), max(ylim)-5, pos=4, font=3, labels="multires spline")

abline(v=60, col="red", lty=2)

plot(lats300, rmse(a.bg, a.aster), type="l", xlab="Latitude",

    ylab="RMSE", ylim=ylim)

lines(lats150, rmse(crop(a.bg, extent(a.srtm)), a.srtm), col="blue")

legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),

    lty=c(1, 1), bty="n")

text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend")

abline(v=60, col="red", lty=2)

# ...with respect to CDEM

plot(lats300, rmse(a.uncor, a.can), type="l", xlab="Latitude",

    ylab="RMSE", ylim=ylim)

text(min(lats300), max(ylim)-5, pos=4, font=3, labels="simple fused")

abline(v=60, col="red", lty=2)

mtext(expression(paste(

    "Aspect discrepancies with respect to Canada DEM (",

    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)

plot(lats300, rmse(a.enblend, a.can), type="l", xlab="Latitude",

    ylab="RMSE", ylim=ylim)

text(min(lats300), max(ylim)-5, pos=4, font=3, labels="multires spline")

abline(v=60, col="red", lty=2)

plot(lats300, rmse(a.bg, a.can), type="l", xlab="Latitude",

    ylab="RMSE", ylim=ylim)

text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend")

abline(v=60, col="red", lty=2)

#

# plot latitudinal profiles of correlation coefficients

#

# simple helper function to calculate row-wise *circular* correlation

# coefficients

corByLat <- function(r1, r2, rows) {

    if (missing(rows)) {

        rows <- 1:nrow(r1)

    }

    m1 <- circular(as.matrix(r1), units="degrees", rotation="clock")

    m2 <- circular(as.matrix(r2), units="degrees", rotation="clock")

    sapply(rows, function(row) {

        p <- cor.circular(m1[row,], m2[row,])

        if (is.null(p)) NA else p

        })

}

par(mfrow=c(2,3), omi=c(1,1,1,1))

ylim <- c(-1, 1)

# ...with respect to ASTER

plot(lats300, corByLat(a.uncor, a.aster), type="l", xlab="Latitude",

    ylab="Circular correlation", ylim=ylim)

lines(lats150, corByLat(crop(a.uncor, extent(a.srtm)), a.srtm), col="blue")

legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),

    lty=c(1, 1), bty="n")

text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused")

abline(v=60, col="red", lty=2)

mtext(expression(paste(

    "Aspect correlations with respect to separate ASTER/SRTM components (",

    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)

plot(lats300, corByLat(a.enblend, a.aster), type="l", xlab="Latitude",

    ylab="Circular correlation", ylim=ylim)

lines(lats150, corByLat(crop(a.enblend, extent(a.srtm)), a.srtm), col="blue")

legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),

    lty=c(1, 1), bty="n")

text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")

abline(v=60, col="red", lty=2)

plot(lats300, corByLat(a.bg, a.aster), type="l", xlab="Latitude",

    ylab="Circular correlation", ylim=ylim)

lines(lats150, corByLat(crop(a.bg, extent(a.srtm)), a.srtm), col="blue")

legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),

    lty=c(1, 1), bty="n")

text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend")

abline(v=60, col="red", lty=2)

# ...with respect to CDEM

plot(lats300, corByLat(a.uncor, a.can), type="l", xlab="Latitude",

    ylab="Circular correlation", ylim=ylim)

text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused")

abline(v=60, col="red", lty=2)

mtext(expression(paste(

    "Aspect correlations with respect to Canada DEM (",

    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)

plot(lats300, corByLat(a.enblend, a.can), type="l", xlab="Latitude",

    ylab="Circular correlation", ylim=ylim)

text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")

abline(v=60, col="red", lty=2)

plot(lats300, corByLat(a.bg, a.can), type="l", xlab="Latitude",

    ylab="Circular correlation", ylim=ylim)

text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend")

abline(v=60, col="red", lty=2)

# close pdf device driver

dev.off()

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

# Quick script to determine whether the filename of each downloaded

# ASTER tile is faithful to its spatial location as specified in the

# GeoTIFF itself.

#

# This script is written to run on eos, based on the file system

# location of downloaded ASTER tiles as of 04-Aug-2011

#

# Jim Regetz

# NCEAS

# Created on 04-Aug-2011

library(rgdal)

# for a given file, return TRUE if file name matches the internally

# specified lower left corner, otherwise return string with those

# coordinates

check <- function(tilename, path=".", silent=TRUE) {

    # build expected filename

    origin <- round(GDALinfo(path.expand(file.path(path, tilename)),

        silent=silent)[c("ll.x", "ll.y")])

    ly <- origin["ll.y"]

    y <- sprintf("%s%02d", if (ly>=0) "N" else "S", abs(ly))

    lx <- origin["ll.x"]

    x <- sprintf("%s%03d", if (lx>=0) "E" else "W", abs(lx))

    expected.name <- paste("ASTGTM_", y, x, "_dem.tif", sep="")

    # compare to actual filename

    if (tilename==expected.name) {

        TRUE

    } else {

        paste(y, x)

    }

}

# produce vector of check results, with the actual file names as vector

# element names (takes ~12 minutes on eos)

aster.dir <- "~organisms/DEM/asterGdem"

aster.tiles <- list.files(aster.dir, pattern="^ASTGTM.*_dem.tif$")

tilecheck <- sapply(aster.tiles, check, path=aster.dir)

# report mismatches

data.frame(expected=tilecheck[tilecheck!="TRUE"])

##                       expected

##ASTGTM_N59E069_dem.tif N63 E109

##ASTGTM_N63E113_dem.tif N69 E107

##ASTGTM_N63E117_dem.tif N69 E113

##ASTGTM_N64E098_dem.tif N70 E117

##ASTGTM_N65E104_dem.tif N73 E084

##ASTGTM_N65E111_dem.tif N73 E098

##ASTGTM_N65E117_dem.tif N66 E130

# GDAL commands for assembling SRTM and ASTER DEM tiles associated with

# the 60N Canada boundary analysis, and resampling into the desired 3"

# (~90m) resolution even-boundary grid.

#

# These commands generate strips of elevation data (GeoTIFFs) along a

# 40-degree longitudinal extent matching (at least one of) Rick Reeves'

# mosaicked CDEM grids. The strips extend 150 pixels south of 60N and

# (in case of ASTER only) 150 pixels north of 60N, which should provide

# a sufficient but not excessive latitudinal range for fixing and

# assessing boundary artifacts.

#

# Note:

#   Working with the original ASTERs yields this warning from GDAL:

#     Warning 1: TIFFReadDirectoryCheckOrder:Invalid TIFF directory;

#     tags are not sorted in ascending order

#   I then ran gdal_translate on several of the offending ASTERs, then

#   repeated the vrt/warp on those -- now without warnings. However, the

#   output data was the same as when I operated on the original files

#   (with warnings), so for the moment I'm just going to ignore the

#   warnings.

#

# Jim Regetz

# NCEAS

# at the time of script creation, these paths were correct on vulcan

export ASTDIR="/home/reeves/active_work/EandO/asterGdem"

export SRTMDIR="/home/reeves/active_work/EandO/CgiarSrtm/SRTM_90m_ASCII_4_1"

# SRTM (also convert to 16bit integer)

gdalbuildvrt srtm.vrt $SRTMDIR/srtm_*_01.asc

gdalwarp -ot Int16 -te -136 59.875 -96 60 -ts 48000 150 -r bilinear \

  srtm.vrt srtm_150below.tif

# ASTER

gdalbuildvrt aster.vrt $ASTDIR/ASTGTM_N59*W*_dem.tif \

  $ASTDIR/ASTGTM_N60*W*_dem.tif

gdalwarp -te -136 59.875 -96 60 -ts 48000 150 -r bilinear \

  aster.vrt aster_150below.tif

gdalwarp -te -136 60 -96 60.125 -ts 48000 150 -r bilinear \

  aster.vrt aster_150above.tif

# note that the top 150 rows of this one are, somewhat surprisingly,

# slightly different from the above!

# gdalwarp -te -136 59.875 -96 60.125 -ts 48000 300 -r bilinear \

#   aster.vrt aster_300straddle.tif

#

# and this yields an even different set of values

# gdalbuildvrt aster_N60.vrt $ASTDIR/ASTGTM_N60*W*_dem.tif

# gdalwarp -te -136 60 -96 60.125 -ts 48000 150 -r bilinear \

#   aster_N60.vrt aster_150above.tif

# R script to apply several kinds of boundary corrections to ASTER/SRTM

# elevation values near the 60N boundary in Canada, and write out new

# GeoTIFFs.

#

# Jim Regetz

# NCEAS

library(raster)

# load relevant SRTM and ASTER data

srtm.south <- raster("srtm_150below.tif")

aster.south <- raster("aster_150below.tif")

aster.north <- raster("aster_150above.tif")

# create difference raster for area of overlap

delta.south <- srtm.south - aster.south

#

# OPTION 1

#

# smooth to the north, by calculating the deltas _at_ the boundary,

# ramping them down to zero with increasing distance from the border,

# and adding them to the north ASTER values

# create simple grid indicating distance (in units of pixels) north from

# boundary, starting at 1 (this is used for both option 1 and option 2)

aster.north.matrix <- as.matrix(aster.north)

ydistN <- nrow(aster.north.matrix) + 1 - row(aster.north.matrix)

# 1a. linear ramp north from SRTM edge

# -- Rick has done this --

# 1b. exponential ramp north from SRTM edge

r <- -0.045

w <- exp(ydistN*r)

aster.north.smooth <- aster.north

aster.north.smooth[] <- values(aster.north) + as.integer(round(t(w) *

    as.matrix(delta.south)[1,]))

writeRaster(aster.north.smooth, file="aster_150above_rampexp.tif")

#

# OPTION 2

#

# smooth to the north, by first using LOESS with values south of 60N to

# model deltas as a function of observed ASTER, then applying the model

# to predict pixel-wise deltas north of 60N, then ramping these

# predicted deltas to zero with increasing distance from the border, and

# adding them to the associated ASTER values

# first fit LOESS on a random subsample of data

# note: doing all the data takes too long, and even doing 50k points

#       seems to be too much for calculating SEs during predict step

set.seed(99)

samp <- sample(ncell(aster.south), 10000)

sampdata <- data.frame(delta=values(delta.south)[samp],

    aster=values(aster.south)[samp])

lo.byaster <- loess(delta ~ aster, data=sampdata)

# now create ASTER prediction grid north of 60N

# TODO: deal with NAs in data (or make sure they are passed through

#       properly in the absence of explicit treatment)?

aster.north.pdelta <- aster.north

aster.north.pdelta[] <- predict(lo.byaster, values(aster.north))

# for actual north ASTER values that exceed the max value used to fit

# LOESS, just use the prediction associated with the maximum

aster.north.pdelta[aster.north<min(sampdata$aster)] <- predict(lo.byaster,

    data.frame(aster=min(sampdata$aster)))

# for actual north ASTER value less than the min value used to fit

# LOESS, just use the prediction associated with the minimum

aster.north.pdelta[aster.north>max(sampdata$aster)] <- predict(lo.byaster,

    data.frame(aster=max(sampdata$aster)))

# 2a: exponential distance-weighting of LOESS predicted deltas

r <- -0.045

w <- exp(r*ydistN)

aster.north.smooth <- aster.north

aster.north.smooth[] <- values(aster.north) + as.integer(round(t(w *

    as.matrix(aster.north.pdelta))))

writeRaster(aster.north.smooth, file="aster_150above_predexp.tif")

# 2b: gaussian distance-weighting of LOESS predicted deltas

r <- -0.001  # weight drops to 0.5 at ~26 cells, ie 2.4km at 3" res

w <- exp(r*ydistN^2)

aster.north.smooth <- aster.north

aster.north.smooth[] <- values(aster.north) + as.integer(round(t(w *

    as.matrix(aster.north.pdelta))))

writeRaster(aster.north.smooth, file="aster_150above_predgau.tif")

#

# OPTION 3

#

# smooth to the south, now by simply taking pixel-wise averages of the

# observed SRTM and ASTER using a distance-based weighting function such

# that the relative contribution of ASTER decays to zero over a few km

# create simple grid indicating distance (in units of pixels) south from

# boundary, starting at 1

aster.south.matrix <- as.matrix(aster.south)

ydistS <- row(aster.south.matrix)

# 3a: gaussian weighting function

r <- -0.001  # weight drops to 0.5 at ~26 cells, or 2.4km at 3" res

w <- exp(-0.001*ydistS^2)

aster.south.smooth <- aster.south

aster.south.smooth[] <- values(srtm.south) - as.integer(round(t(w *

    as.matrix(delta.south))))

aster.south.smooth[aster.south.smooth<0] <- 0

writeRaster(aster.south.smooth, file="dem_150below_blendgau.tif")

# GDAL commands to produced fused DEMs in the vicinity of the 60N Canada

# boundary, using several "boundary-corrected" variants as well as the

# original uncorrected DEMs. Note that the multiresolution spline is not

# included here, because the associated fused layer is already produced

# in its entirety during that process.

#

# Jim Regetz

# NCEAS

# uncorrected fused layer

gdalwarp -ot Int16 -te -136 59.875 -96 60.125 -ts 48000 300 \

  srtm_150below.tif aster_150above.tif fused_300straddle.tif

# exponential ramp of boundary delta to the north

gdalwarp -ot Int16 -te -136 59.875 -96 60.125 -ts 48000 300 \

  srtm_150below.tif aster_150above_rampexp.tif fused_300straddle_rampexp.tif

# exponential blend of predicted deltas to the north

gdalwarp -ot Int16 -te -136 59.875 -96 60.125 -ts 48000 300 \

  srtm_150below.tif aster_150above_predexp.tif fused_300straddle_predexp.tif

# gaussian blend of predicted deltas to the north

gdalwarp -ot Int16 -te -136 59.875 -96 60.125 -ts 48000 300 \

  srtm_150below.tif aster_150above_predgau.tif fused_300straddle_predgau.tif

# gaussian blend of SRTM/ASTER to the south

gdalwarp -ot Int16 -te -136 59.875 -96 60.125 -ts 48000 300 \

  dem_150below_blendgau.tif aster_150above.tif fused_300straddle_blendgau.tif

# GDAL commands to generate hillshade images using the different fused

# DEMs we've produced in the 60N Canada boundary region.

#

# Also included below is a bit of R code to generate a PNG of quick

# hillshade visuals of two 1-degree wide swaths (one in the more

# mountainous west, and one in the flatter east) spanning the 60N

# boundary in Canada, for several different fused layers.

#

# Jim Regetz

# NCEAS

export DEMDIR=~/media/temp/terrain/dem

export HSDIR=~/media/temp/terrain/hillshade

gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle.tif \

  $HSDIR/fused_300straddle_h.tif

gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_rampexp.tif \

  $HSDIR/fused_300straddle_rampexp_h.tif

gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_predexp.tif \

  $HSDIR/fused_300straddle_predexp_h.tif

gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_predgau.tif \

  $HSDIR/fused_300straddle_predgau_h.tif

gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_blendgau.tif \

  $HSDIR/fused_300straddle_blendgau_h.tif

gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_enblend.tif \

  $HSDIR/fused_300straddle_enblend_h.tif

# Use R to generate a PDF with some visuals 

echo '

  library(raster)

  hsdir <- "/home/regetz/media/temp/terrain/hillshade"

  h.uncor <- raster(file.path(hsdir, "fused_300straddle_h.tif"))

  h.re <- raster(file.path(hsdir, "fused_300straddle_rampexp_h.tif"))

  h.enblend <- raster(file.path(hsdir, "fused_300straddle_enblend_h.tif"))

  h.bg <- raster(file.path(hsdir, "fused_300straddle_blendgau_h.tif"))

  window1 <- extent(-135, -134, 59.875, 60.125)

  window2 <- extent(-118, -117, 59.875, 60.125)

  png("boundary-hillshade.png", height=10, width=7.5, units="in", res=600)

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

  plot(crop(h.uncor, window1), legend=FALSE)

  mtext("Simple fuse (hillshade)", adj=0, cex=0.8)

  plot(crop(h.uncor, window2), legend=FALSE)

  plot(crop(h.re, window1), legend=FALSE)

  mtext("North exponential ramp (hillshade)", adj=0, cex=0.8)

  plot(crop(h.re, window2), legend=FALSE)

  plot(crop(h.enblend, window1), legend=FALSE)

  mtext("Multiresolution spline (hillshade)", adj=0, cex=0.8)

  plot(crop(h.enblend, window2), legend=FALSE)

  plot(crop(h.bg, window1), legend=FALSE)

  mtext("Gaussian weighted average (hillshade)", adj=0, cex=0.8)

  plot(crop(h.bg, window2), legend=FALSE)

  dev.off()

' | Rscript --vanilla -

# R code to plot latitudinal profiles of mean elevation, along with both

# RMSE and correlation coefficients comparing fused layers with both the

# raw ASTER and with the Canada DEM

#

# Jim Regetz

# NCEAS

# Created on 08-Jun-2011

library(raster)

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

# load elevation rasters

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

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

d.uncor <- raster(file.path(datadir, "fused_300straddle.tif"))

d.enblend <- raster(file.path(datadir, "fused_300straddle_enblend.tif"))

d.bg <- raster(file.path(datadir, "fused_300straddle_blendgau.tif"))

d.can <- raster(file.path(datadir, "cdem_300straddle.tif"))

# extract raster latitudes for later

lats300 <- yFromRow(d.aster, 1:nrow(d.aster))

lats150 <- yFromRow(d.srtm, 1:nrow(d.srtm))

#

# plot latitudinal profiles of mean elevation

#

# initialize output pdf device driver

pdf("elevation-assessment.pdf", height=8, width=11.5)

par(mfrow=c(2,2), omi=c(1,1,1,1))

ylim <- c(540, 575)

plot(lats300, rowMeans(as.matrix(d.can), na.rm=TRUE), type="l",

    xlab="Latitude", ylab="Mean elevation", ylim=ylim)

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Original DEMs")

lines(lats300, rowMeans(as.matrix(d.aster), na.rm=TRUE), col="blue")

lines(lats150, rowMeans(as.matrix(d.srtm), na.rm=TRUE), col="red")

legend("bottomright", legend=c("SRTM", "CDED", "ASTER"), col=c("red",

    "black", "blue"), lty=c(1, 1), bty="n")

abline(v=60, col="red", lty=2)

mtext(expression(paste("Latitudinal profiles of mean elevation (",

    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)

plot(lats300, rowMeans(as.matrix(d.uncor), na.rm=TRUE), type="l",

    xlab="Latitude", ylab="Mean elevation", ylim=ylim)

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="simple fuse")

abline(v=60, col="red", lty=2)

plot(lats300, rowMeans(as.matrix(d.enblend), na.rm=TRUE), type="l",

    xlab="Latitude", ylab="Mean elevation", ylim=ylim)

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="multires spline")

abline(v=60, col="red", lty=2)

plot(lats300, rowMeans(as.matrix(d.bg), na.rm=TRUE), type="l",

    xlab="Latitude", ylab="Mean elevation", ylim=ylim)

text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend")

abline(v=60, col="red", lty=2)

#

# plot latitudinal profiles of RMSE

#

# simple helper function to calculate row-wise RMSEs