# 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

# Created on 08-Jun-2011

# 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

# Created on 08-Jun-2011

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

#

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