1 |
3afd1c7b
|
Jim Regetz
|
# 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 |
cc5360b4
|
Jim Regetz
|
# 1a. linear ramp north from SRTM edge
|
32 |
|
|
# -- Rick has done this --
|
33 |
3afd1c7b
|
Jim Regetz
|
|
34 |
cc5360b4
|
Jim Regetz
|
# 1b. exponential ramp north from SRTM edge
|
35 |
3afd1c7b
|
Jim Regetz
|
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")
|