|
1
|
# R code to smooth the boundary area N59 to N60 to the south, by simply taking pixel-wise
|
|
2
|
# averages of the observed SRTM and ASTER using a distance-based weighting function such
|
|
3
|
# that the relative contribution of ASTER decays to zero over a few km
|
|
4
|
# create simple grid indicating distance (in units of pixels) south from
|
|
5
|
# boundary , starting at 1.
|
|
6
|
#
|
|
7
|
#
|
|
8
|
# Original author: Jim Regetz
|
|
9
|
# [13-JUL-2011]
|
|
10
|
# Edits by Yuni to apply ASTER GDEM2, instead of GDEM1, as well as auto-applying
|
|
11
|
# for each nine region. [9-Dec-2011]
|
|
12
|
#
|
|
13
|
|
|
14
|
|
|
15
|
|
|
16
|
|
|
17
|
library(raster)
|
|
18
|
|
|
19
|
aster.dir <- "/data/project/organisms/DEM/Yuni/Data/aster2"
|
|
20
|
srtm.dir <- "/data/project/organisms/DEM/Yuni/Data/srtm"
|
|
21
|
aster.below <- list.files(aster.dir, pattern="^aster2_.*_below.tif")
|
|
22
|
aster.above <- list.files(aster.dir, pattern="^aster2_.*_above.tif")
|
|
23
|
srtm.below <- list.files(srtm.dir, pattern="^srtm_.*_below.tif$")
|
|
24
|
|
|
25
|
|
|
26
|
blendgau <- function(aster.tile.below, aster.tile.above, srtm.tile.below, aster.dir, srtm.dir){
|
|
27
|
# load relevant SRTM and ASTER data
|
|
28
|
srtm.south <- raster(file.path(srtm.dir, srtm.tile.below))
|
|
29
|
aster.south <- raster(file.path(aster.dir, aster.tile.below))
|
|
30
|
aster.north <- raster(file.path(aster.dir, aster.tile.above))
|
|
31
|
|
|
32
|
# create difference raster for area of overlap
|
|
33
|
delta.south <- srtm.south - aster.south
|
|
34
|
|
|
35
|
aster.south.matrix <- as.matrix(aster.south)
|
|
36
|
ydistS <- row(aster.south.matrix)
|
|
37
|
|
|
38
|
r <- -0.001 # weight drops to 0.5 at ~26 cells, or 2.4km at 3" res
|
|
39
|
w <- exp(-0.001*ydistS^2)
|
|
40
|
aster.south.smooth <- aster.south
|
|
41
|
aster.south.smooth[] <- values(srtm.south) - as.integer(round(t(w *
|
|
42
|
as.matrix(delta.south))))
|
|
43
|
aster.south.smooth[aster.south.smooth <0] <- 0
|
|
44
|
file.name <- paste("aster2_", substr(srtm.tile.below, 6, 13), "_below_blendgau.tif", sep="")
|
|
45
|
writeRaster(aster.south.smooth, file.path(aster.dir, file.name))
|
|
46
|
}
|
|
47
|
|
|
48
|
|
|
49
|
|
|
50
|
|
|
51
|
num.regions <- 9
|
|
52
|
makeGau <-lapply(1:num.regions, function(region, aster.dir, srtm.dir){
|
|
53
|
aster.tile.below <- aster.below[region]
|
|
54
|
aster.tile.above <- aster.above[region]
|
|
55
|
srtm.tile.below <- srtm.below[region]
|
|
56
|
blendgau(aster.tile.below, aster.tile.above, srtm.tile.below, aster.dir, srtm.dir)
|
|
57
|
}, aster.dir=aster.dir, srtm.dir=srtm.dir)
|
|
58
|
|
|
59
|
|
|
60
|
|
