1
|
# R script to apply several kinds of boundary corrections to ASTER/SRTM
|
2
|
# elevation values near the 60N boundary in Canada, and write out new
|
3
|
# GeoTIFFs.
|
4
|
#
|
5
|
# Jim Regetz
|
6
|
# NCEAS
|
7
|
|
8
|
library(raster)
|
9
|
|
10
|
# load relevant SRTM and ASTER data
|
11
|
srtm.south <- raster("srtm_150below.tif")
|
12
|
aster.south <- raster("aster_150below.tif")
|
13
|
aster.north <- raster("aster_150above.tif")
|
14
|
|
15
|
# create difference raster for area of overlap
|
16
|
delta.south <- srtm.south - aster.south
|
17
|
|
18
|
#
|
19
|
# OPTION 1
|
20
|
#
|
21
|
|
22
|
# smooth to the north, by calculating the deltas _at_ the boundary,
|
23
|
# ramping them down to zero with increasing distance from the border,
|
24
|
# and adding them to the north ASTER values
|
25
|
|
26
|
# create simple grid indicating distance (in units of pixels) north from
|
27
|
# boundary, starting at 1 (this is used for both option 1 and option 2)
|
28
|
aster.north.matrix <- as.matrix(aster.north)
|
29
|
ydistN <- nrow(aster.north.matrix) + 1 - row(aster.north.matrix)
|
30
|
|
31
|
# 1a. linear ramp north from SRTM edge
|
32
|
# -- Rick has done this --
|
33
|
|
34
|
# 1b. exponential ramp north from SRTM edge
|
35
|
r <- -0.045
|
36
|
w <- exp(ydistN*r)
|
37
|
aster.north.smooth <- aster.north
|
38
|
aster.north.smooth[] <- values(aster.north) + as.integer(round(t(w) *
|
39
|
as.matrix(delta.south)[1,]))
|
40
|
writeRaster(aster.north.smooth, file="aster_150above_rampexp.tif")
|
41
|
|
42
|
#
|
43
|
# OPTION 2
|
44
|
#
|
45
|
|
46
|
# smooth to the north, by first using LOESS with values south of 60N to
|
47
|
# model deltas as a function of observed ASTER, then applying the model
|
48
|
# to predict pixel-wise deltas north of 60N, then ramping these
|
49
|
# predicted deltas to zero with increasing distance from the border, and
|
50
|
# adding them to the associated ASTER values
|
51
|
|
52
|
# first fit LOESS on a random subsample of data
|
53
|
# note: doing all the data takes too long, and even doing 50k points
|
54
|
# seems to be too much for calculating SEs during predict step
|
55
|
set.seed(99)
|
56
|
samp <- sample(ncell(aster.south), 10000)
|
57
|
sampdata <- data.frame(delta=values(delta.south)[samp],
|
58
|
aster=values(aster.south)[samp])
|
59
|
lo.byaster <- loess(delta ~ aster, data=sampdata)
|
60
|
|
61
|
# now create ASTER prediction grid north of 60N
|
62
|
# TODO: deal with NAs in data (or make sure they are passed through
|
63
|
# properly in the absence of explicit treatment)?
|
64
|
aster.north.pdelta <- aster.north
|
65
|
aster.north.pdelta[] <- predict(lo.byaster, values(aster.north))
|
66
|
# for actual north ASTER values that exceed the max value used to fit
|
67
|
# LOESS, just use the prediction associated with the maximum
|
68
|
aster.north.pdelta[aster.north<min(sampdata$aster)] <- predict(lo.byaster,
|
69
|
data.frame(aster=min(sampdata$aster)))
|
70
|
# for actual north ASTER value less than the min value used to fit
|
71
|
# LOESS, just use the prediction associated with the minimum
|
72
|
aster.north.pdelta[aster.north>max(sampdata$aster)] <- predict(lo.byaster,
|
73
|
data.frame(aster=max(sampdata$aster)))
|
74
|
|
75
|
# 2a: exponential distance-weighting of LOESS predicted deltas
|
76
|
r <- -0.045
|
77
|
w <- exp(r*ydistN)
|
78
|
aster.north.smooth <- aster.north
|
79
|
aster.north.smooth[] <- values(aster.north) + as.integer(round(t(w *
|
80
|
as.matrix(aster.north.pdelta))))
|
81
|
writeRaster(aster.north.smooth, file="aster_150above_predexp.tif")
|
82
|
|
83
|
# 2b: gaussian distance-weighting of LOESS predicted deltas
|
84
|
r <- -0.001 # weight drops to 0.5 at ~26 cells, ie 2.4km at 3" res
|
85
|
w <- exp(r*ydistN^2)
|
86
|
aster.north.smooth <- aster.north
|
87
|
aster.north.smooth[] <- values(aster.north) + as.integer(round(t(w *
|
88
|
as.matrix(aster.north.pdelta))))
|
89
|
writeRaster(aster.north.smooth, file="aster_150above_predgau.tif")
|
90
|
|
91
|
#
|
92
|
# OPTION 3
|
93
|
#
|
94
|
|
95
|
# smooth to the south, now by simply taking pixel-wise averages of the
|
96
|
# observed SRTM and ASTER using a distance-based weighting function such
|
97
|
# that the relative contribution of ASTER decays to zero over a few km
|
98
|
|
99
|
# create simple grid indicating distance (in units of pixels) south from
|
100
|
# boundary, starting at 1
|
101
|
aster.south.matrix <- as.matrix(aster.south)
|
102
|
ydistS <- row(aster.south.matrix)
|
103
|
|
104
|
# 3a: gaussian weighting function
|
105
|
r <- -0.001 # weight drops to 0.5 at ~26 cells, or 2.4km at 3" res
|
106
|
w <- exp(-0.001*ydistS^2)
|
107
|
aster.south.smooth <- aster.south
|
108
|
aster.south.smooth[] <- values(srtm.south) - as.integer(round(t(w *
|
109
|
as.matrix(delta.south))))
|
110
|
aster.south.smooth[aster.south.smooth<0] <- 0
|
111
|
writeRaster(aster.south.smooth, file="dem_150below_blendgau.tif")
|