# Snippets of GDAL commands and R code for processing DEMs

# Jim Regetz

# Created on 08-Jun-2011

#

# 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 first ran gdal_translate on each of the ASTERs, then repeated the

# vrt/warp on those (without warnings), but the output 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?

#=======================================================================

# bash -- resample source DEMs into desired grids near the 60N boundary

#=======================================================================

# generate strips of data along a 40-degree longitudinal extent matching

# (at least one of) Rick's mosaicked CDEM grids; strips extend 150

# pixels south of border and (in case of aster) north of border

# these are currently 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 -- apply several kinds of boundary fixes and write out new GeoTIFFs

#=======================================================================

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)

# 1b. linear ramp north from SRTM edge

# -- Rick is doing this --

# 2b. exponential ramp north from SRTM edge

# -- Rick is also doing this, but here it is... --

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(ydistN*r)

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(-0.001*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")

#=======================================================================

# bash -- fuse DEMS, generate hillshade

#=======================================================================

#

# create simple fused layers

#

# 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

#

# hillshade the different fused DEMs created above

#

gdaldem hillshade -s 111120 fused_300straddle.tif fused_300straddle_hs.tif

gdaldem hillshade -s 111120 fused_300straddle_rampexp.tif fused_300straddle_rampexp_hs.tif

gdaldem hillshade -s 111120 fused_300straddle_predexp.tif fused_300straddle_predexp_hs.tif

gdaldem hillshade -s 111120 fused_300straddle_predgau.tif fused_300straddle_predgau_hs.tif

gdaldem hillshade -s 111120 fused_300straddle_blendgau.tif fused_300straddle_blendgau_hs.tif

#=======================================================================

# R -- generate some quick hillshade visuals of a 1-degree wide swath

#=======================================================================

library(raster)

uncorrected <- raster("fused_300straddle_hs.tif")

rampexp <- raster("fused_300straddle_rampexp_hs.tif")

blendgau <- raster("fused_300straddle_blendgau_hs.tif")

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

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

par(mfrow=c(3,1))

plot(crop(uncorrected, window), main="uncorrected (hillshade)")

plot(crop(rampexp, window), main="north exponential ramp (hillshade)")

plot(crop(blendgau, window), main="south gaussian blend (hillshade)")

dev.off()

#=======================================================================

# R -- assess boundary artifacts with respect to slope

#=======================================================================

s.aster <- raster("aster_300straddle_s.tif")

s.srtm <- raster("srtm_150below_s.tif")

s.uncor <- raster("fused_300straddle_s.tif")

s.eramp <- raster("fused_300straddle_rampexp_s.tif")

s.bg <- raster("fused_300straddle_blendgau_s.tif")

s.can <- raster("cdem_300straddle_s.tif")

rmse <- function(r1, r2) {

    sqrt(rowMeans(as.matrix((r1 - r2)^2), na.rm=TRUE))

}

pdf("slope-rmse.pdf", height=8, width=11.5)

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

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

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

# Latitudinal RMSE profiles with respect to ASTER

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

    ylab="RMSE", ylim=c(0, 5))

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

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

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

text(min(lats300), 4.5, pos=4, font=3, labels="uncorrected")

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

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

    adj=0, line=2, font=2)

plot(lats300, rmse(s.eramp, s.aster), type="l", xlab="Latitude",

    ylab="RMSE", ylim=c(0, 5))

lines(lats150, rmse(crop(s.eramp, extent(s.srtm)), s.srtm), col="blue")

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

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

text(min(lats300), 4.5, pos=4, font=3, labels="exponential ramp")

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

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

    ylab="RMSE", ylim=c(0, 5))

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

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

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

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

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

# Latitudinal RMSE profiles with respect to CDEM

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

    ylab="RMSE", ylim=c(0, 5))

text(min(lats300), 4.5, pos=4, font=3, labels="uncorrected")

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

mtext("Slope discrepancies with respect to Canada DEM",

    adj=0, line=2, font=2)

plot(lats300, rmse(s.eramp, s.can), type="l", xlab="Latitude",

    ylab="RMSE", ylim=c(0, 5))

text(min(lats300), 4.5, pos=4, font=3, labels="exponential ramp")

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

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

    ylab="RMSE", ylim=c(0, 5))

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

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

dev.off()

####

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

    if (missing(rows)) {

        rows <- 1:nrow(r1)

    }

    m1 <- as.matrix(r1)

    m2 <- as.matrix(r2)

    sapply(rows, function(row) cor(m1[row,], m2[row,],

        use="pairwise.complete.obs"))

}

pdf("slope-corr.pdf", height=8, width=11.5)

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

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

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

ylim <- c(0.65, 1)

# Latitudinal RMSE profiles with respect to ASTER

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

    ylab="Correlation", ylim=ylim)

lines(lats150, corByLat(crop(s.uncor, extent(s.srtm)), s.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="uncorrected")

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

mtext("Slope correlations with separate ASTER/SRTM components",

    adj=0, line=2, font=2)

plot(lats300, corByLat(s.eramp, s.aster), type="l", xlab="Latitude",

    ylab="Correlation", ylim=ylim)

lines(lats150, corByLat(crop(s.eramp, extent(s.srtm)), s.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="exponential ramp")

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

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

    ylab="Correlation", ylim=ylim)

lines(lats150, corByLat(crop(s.bg, extent(s.srtm)), s.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)

# Latitudinal correlation profiles with respect to CDEM

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

    ylab="Correlation", ylim=ylim)

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

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

mtext("Slope correlations with Canada DEM",

    adj=0, line=2, font=2)

plot(lats300, corByLat(s.eramp, s.can), type="l", xlab="Latitude",

    ylab="Correlation", ylim=ylim)

text(min(lats300), min(ylim), pos=4, font=3, labels="exponential ramp")

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

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

    ylab="Correlation", ylim=ylim)

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

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

dev.off()