|
1
|
# Modificaion of 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 (organisms/DEM/Yuni/scripts/toProduceData/gaussian.r)
|
|
9
|
# [13-JUL-2011]
|
|
10
|
# Edits by Natalie to apply to resampled ASTER GDEM2, in 5 x 5 grids with no 1/2 pixel offset, and
|
|
11
|
# shifterd SRTM tiles, with no 1/2 pixel offset. [2-Feb-2012]
|
|
12
|
# NOTE- A similar version of this script exists in "organisms/DEM/Yuni/scripts/toProduceData/gaussian.r" as
|
|
13
|
# written and edited by Jim and Yuni.
|
|
14
|
|
|
15
|
|
|
16
|
library(raster)
|
|
17
|
|
|
18
|
#Set working directories and filenames
|
|
19
|
aster.dir <- "/data/project/organisms/DEM/asterGdem2/90m_NoPixelOffset/Mosaiced/N59to60"
|
|
20
|
srtm.dir <- "/data/project/organisms/DEM/cgiarSrtm/SRTM_90m_ASCII_4_1/Tiles_Resampled/Mosaiced"
|
|
21
|
out.dir <- "/data/project/organisms/DEM/GlobalProduct"
|
|
22
|
aster.below.East <- paste (aster.dir, "/N59to60E000_180.tif", sep="")
|
|
23
|
srtm.East <- paste (srtm.dir, "/SRTM_N59to60E000_180.tif",sep="")
|
|
24
|
aster.below.West <- paste(aster.dir, "/N59to60W180_000.tif",sep="")
|
|
25
|
srtm.West <- paste(srtm.dir, "/SRTM_N59to60W180_000.tif",sep="")
|
|
26
|
|
|
27
|
blendgau <- function(aster.tile, srtm.tile, aster.dir, srtm.dir, out.dir, hemi){
|
|
28
|
# load relevant SRTM and ASTER data
|
|
29
|
srtm.south <- raster(srtm.tile)
|
|
30
|
aster.south <- raster(aster.tile)
|
|
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("/Blended_N59to60_", hemi, ".tif", sep="")
|
|
45
|
writeRaster(aster.south.smooth, file.path(out.dir, file.name))
|
|
46
|
}
|
|
47
|
|
|
48
|
blendgau(aster.below.East, srtm.East, aster.dir, srtm.dir, out.dir, hemi="East")
|
|
49
|
blendgau(aster.below.West, srtm.West, aster.dir, srtm.dir, out.dir, hemi="West")
|
