| 1 |
f12b574f
|
Jim Regetz
|
# Code to produce enblended DEM (i.e., using multiresolution splines as
|
| 2 |
933632ec
|
Jim Regetz
|
# described by Burt & Adelson 1983) in the 60N boundary region. After
|
| 3 |
|
|
# appropriately preparing the SRTM and ASTER layers, this code makes a
|
| 4 |
|
|
# system call to run 'enblend' (v. 4.0) on the inputs, then
|
| 5 |
|
|
# post-processes the resulting image to yield a single band geotiff with
|
| 6 |
|
|
# datatype of 16bit signed integer, matching the input data.
|
| 7 |
f12b574f
|
Jim Regetz
|
#
|
| 8 |
933632ec
|
Jim Regetz
|
# Somewhat arbitrarily, in the code below the input DEMs are prepped
|
| 9 |
|
|
# such that the area of SRTM/ASTER overlap is the first 75 rows below
|
| 10 |
|
|
# 60N (i.e., a zone extending ~6.75km south of the boundary).
|
| 11 |
f12b574f
|
Jim Regetz
|
#
|
| 12 |
|
|
# Jim Regetz
|
| 13 |
|
|
# NCEAS
|
| 14 |
|
|
# Created on 29-Jun-2011
|
| 15 |
|
|
|
| 16 |
|
|
library(raster)
|
| 17 |
|
|
|
| 18 |
|
|
demdir <- "/home/regetz/media/temp/terrain/dem"
|
| 19 |
|
|
|
| 20 |
|
|
# read in aster data
|
| 21 |
|
|
aster <- raster(file.path(demdir, "aster_300straddle.tif"))
|
| 22 |
|
|
# create alpha layer for aster
|
| 23 |
|
|
alpha <- aster
|
| 24 |
|
|
alpha[] <- 0
|
| 25 |
|
|
alpha[1:225, ] <- 1
|
| 26 |
|
|
# write enblend-ready tif out to disk
|
| 27 |
|
|
writeRaster(brick(aster, alpha), file="aster-enblend.tif",
|
| 28 |
|
|
options="ALPHA=YES")
|
| 29 |
|
|
|
| 30 |
|
|
# read in srtm
|
| 31 |
|
|
srtm.lower <- raster("../srtm_150below.tif")
|
| 32 |
|
|
# expand to match full image (i.e., same extent as aster), holding srtm
|
| 33 |
|
|
# data below 60N and nothing above
|
| 34 |
|
|
srtm <- aster
|
| 35 |
|
|
srtm[] <- NA
|
| 36 |
|
|
srtm[151:300, ] <- values(srtm.lower)
|
| 37 |
|
|
# create alpha layer for srtm
|
| 38 |
|
|
alpha <- srtm
|
| 39 |
|
|
alpha[] <- 0
|
| 40 |
|
|
alpha[151:300, ] <- 1
|
| 41 |
|
|
writeRaster(brick(srtm, alpha), file="srtm-enblend.tif",
|
| 42 |
|
|
options="ALPHA=YES")
|
| 43 |
|
|
|
| 44 |
933632ec
|
Jim Regetz
|
# run 'enblend'
|
| 45 |
f12b574f
|
Jim Regetz
|
system(paste("enblend --verbose=6 -o enblend.tif",
|
| 46 |
|
|
"aster-enblend.tif srtm-enblend.tif"))
|
| 47 |
|
|
|
| 48 |
|
|
# post-process enblended DEM
|
| 49 |
|
|
e <- raster("enblend.tif")
|
| 50 |
|
|
e2 <- aster
|
| 51 |
933632ec
|
Jim Regetz
|
# round to nearest integer, and write out the result as a geotiff
|
| 52 |
4bc2a9e9
|
Jim Regetz
|
e2[] <- as.integer(round(values(e), 0))
|
| 53 |
f12b574f
|
Jim Regetz
|
writeRaster(e2, file=file.path(demdir, "fused_300straddle_enblend.tif"),
|
| 54 |
|
|
options="COMPRESS=NONE", datatype="INT2S")
|
