\author{Jim Regetz, NCEAS}

\date{Last update: 13 Jul 2011}

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

 \item SRTM vs ASTER differences

   \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 Fusion via northward exponential rampdown of boundary delta

   \item eliminates elevation cliff at 60\N

   \item introduces north-south ridging artifacts

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

 \item Fusion via multiresolution spline

   \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

 \item Fusion via Gaussian weighted average of SRTM/ASTER

   \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

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

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.

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

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

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

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

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

\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

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

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

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.

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

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)}, 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 STRM (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?

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

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

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

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

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

(muliresolution spline and Gaussian weighted average) exhibit a

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

throughout this zone. 

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.

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

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

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

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

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

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

were calculated for each latitudinal band. Flow direction correlations

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

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

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

  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}

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

  observed median difference between SRTM and ASTER (\Sexpr{delta.median}m).

  Inset histogram shows distribution of differences, excluding absolute

  differences >60.]{

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

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

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

\clearpage

\appendix

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

\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

    flow directions using the various fused elevation rasters and their

    components DEMS as inputs.},

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