Project

General

Profile

Download (10.3 KB) Statistics
| Branch: | Revision:
1 b02e39a0 Jim Regetz
# R code to plot latitudinal profiles of mean elevation, along with both
2
# RMSE and correlation coefficients comparing fused layers with both the
3
# raw ASTER and with the Canada DEM
4
#
5
# Jim Regetz
6
# NCEAS
7
# Created on 08-Jun-2011
8
9
library(raster)
10
11
datadir <- "/home/regetz/media/temp/terrain/dem"
12
13
# load elevation rasters
14
d.aster <- raster(file.path(datadir, "aster_300straddle.tif"))
15
d.srtm <- raster(file.path(datadir, "srtm_150below.tif"))
16
d.uncor <- raster(file.path(datadir, "fused_300straddle.tif"))
17 7cf747b0 Jim Regetz
d.enblend <- raster(file.path(datadir, "fused_300straddle_enblend.tif"))
18 b02e39a0 Jim Regetz
d.bg <- raster(file.path(datadir, "fused_300straddle_blendgau.tif"))
19
d.can <- raster(file.path(datadir, "cdem_300straddle.tif"))
20
21
# extract raster latitudes for later
22
lats300 <- yFromRow(d.aster, 1:nrow(d.aster))
23
lats150 <- yFromRow(d.srtm, 1:nrow(d.srtm))
24
25
26
#
27
# plot latitudinal profiles of mean elevation
28
#
29
30 6626e7db Jim Regetz
# initialize output pdf device driver
31
pdf("elevation-assessment.pdf", height=8, width=11.5)
32
33 b02e39a0 Jim Regetz
par(mfrow=c(2,2), omi=c(1,1,1,1))
34
35 7cf747b0 Jim Regetz
ylim <- c(540, 575)
36 b02e39a0 Jim Regetz
37 7cf747b0 Jim Regetz
plot(lats300, rowMeans(as.matrix(d.can), na.rm=TRUE), type="l",
38 b02e39a0 Jim Regetz
    xlab="Latitude", ylab="Mean elevation", ylim=ylim)
39 7cf747b0 Jim Regetz
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Original DEMs")
40
lines(lats300, rowMeans(as.matrix(d.aster), na.rm=TRUE), col="blue")
41
lines(lats150, rowMeans(as.matrix(d.srtm), na.rm=TRUE), col="red")
42
legend("bottomright", legend=c("SRTM", "CDED", "ASTER"), col=c("red",
43
    "black", "blue"), lty=c(1, 1), bty="n")
44 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
45
mtext(expression(paste("Latitudinal profiles of mean elevation (",
46
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
47
48 7cf747b0 Jim Regetz
plot(lats300, rowMeans(as.matrix(d.uncor), na.rm=TRUE), type="l",
49 b02e39a0 Jim Regetz
    xlab="Latitude", ylab="Mean elevation", ylim=ylim)
50 7cf747b0 Jim Regetz
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="simple fuse")
51 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
52
53 7cf747b0 Jim Regetz
plot(lats300, rowMeans(as.matrix(d.enblend), na.rm=TRUE), type="l",
54 b02e39a0 Jim Regetz
    xlab="Latitude", ylab="Mean elevation", ylim=ylim)
55 7cf747b0 Jim Regetz
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="multires spline")
56 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
57
58
plot(lats300, rowMeans(as.matrix(d.bg), na.rm=TRUE), type="l",
59
    xlab="Latitude", ylab="Mean elevation", ylim=ylim)
60
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend")
61
abline(v=60, col="red", lty=2)
62
63
64
#
65
# plot latitudinal profiles of RMSE
66
#
67
68
# simple helper function to calculate row-wise RMSEs
69
rmse <- function(r1, r2, na.rm=TRUE, use) {
70
    diffs <- abs(as.matrix(r1) - as.matrix(r2))
71
    if (!missing(use)) diffs[!use] <- NA
72
    sqrt(rowMeans(diffs^2, na.rm=na.rm))
73
}
74
75
par(mfrow=c(2,3), omi=c(1,1,1,1))
76
77
ylim <- c(0, 35)
78
79
# ...with respect to ASTER
80
plot(lats300, rmse(d.uncor, d.aster), type="l", xlab="Latitude",
81
    ylab="RMSE", ylim=ylim)
82
lines(lats150, rmse(crop(d.uncor, extent(d.srtm)), d.srtm), col="blue")
83
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
84
    lty=c(1, 1), bty="n")
85 7cf747b0 Jim Regetz
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="simple fuse")
86 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
87
mtext(expression(paste(
88
    "Elevation discrepancies with respect to separate ASTER/SRTM components (",
89
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
90
91 7cf747b0 Jim Regetz
plot(lats300, rmse(d.enblend, d.aster), type="l", xlab="Latitude",
92 b02e39a0 Jim Regetz
    ylab="RMSE", ylim=ylim)
93 7cf747b0 Jim Regetz
lines(lats150, rmse(crop(d.enblend, extent(d.srtm)), d.srtm), col="blue")
94 b02e39a0 Jim Regetz
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
95
    lty=c(1, 1), bty="n")
96 7cf747b0 Jim Regetz
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="multires spline")
97 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
98
99
plot(lats300, rmse(d.bg, d.aster), type="l", xlab="Latitude",
100
    ylab="RMSE", ylim=ylim)
101
lines(lats150, rmse(crop(d.bg, extent(d.srtm)), d.srtm), col="blue")
102
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
103
    lty=c(1, 1), bty="n")
104
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="gaussian blend")
105
abline(v=60, col="red", lty=2)
106
107
# ...with respect to CDEM
108
plot(lats300, rmse(d.uncor, d.can), type="l", xlab="Latitude",
109
    ylab="RMSE", ylim=ylim)
110 7cf747b0 Jim Regetz
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="simple fuse")
111 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
112
mtext(expression(paste(
113
    "Elevation discrepancies with respect to Canada DEM (",
114
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
115
116 7cf747b0 Jim Regetz
plot(lats300, rmse(d.enblend, d.can), type="l", xlab="Latitude",
117 b02e39a0 Jim Regetz
    ylab="RMSE", ylim=ylim)
118 7cf747b0 Jim Regetz
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="multires spline")
119 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
120
121
plot(lats300, rmse(d.bg, d.can), type="l", xlab="Latitude",
122
    ylab="RMSE", ylim=ylim)
123
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="gaussian blend")
124
abline(v=60, col="red", lty=2)
125
126
127
#
128
# plot latitudinal profiles of correlation coefficients
129
#
130
131
# simple helper function to calculate row-wise correlation coefficients
132
corByLat <- function(r1, r2, rows) {
133
    if (missing(rows)) {
134
        rows <- 1:nrow(r1)
135
    }
136
    m1 <- as.matrix(r1)
137
    m2 <- as.matrix(r2)
138
    sapply(rows, function(row) cor(m1[row,], m2[row,],
139
        use="pairwise.complete.obs"))
140
}
141
142
par(mfrow=c(2,3), omi=c(1,1,1,1))
143
144
ylim <- c(0.99, 1)
145
146
# ...with respect to ASTER
147
plot(lats300, corByLat(d.uncor, d.aster), type="l", xlab="Latitude",
148
    ylab="Correlation", ylim=ylim)
149
lines(lats150, corByLat(crop(d.uncor, extent(d.srtm)), d.srtm), col="blue")
150
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
151
    lty=c(1, 1), bty="n")
152 7cf747b0 Jim Regetz
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fuse")
153 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
154
mtext(expression(paste(
155
    "Elevation correlations with respect to separate ASTER/SRTM components (",
156
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
157
158 7cf747b0 Jim Regetz
plot(lats300, corByLat(d.enblend, d.aster), type="l", xlab="Latitude",
159 b02e39a0 Jim Regetz
    ylab="Correlation", ylim=ylim)
160 7cf747b0 Jim Regetz
lines(lats150, corByLat(crop(d.enblend, extent(d.srtm)), d.srtm), col="blue")
161 b02e39a0 Jim Regetz
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
162
    lty=c(1, 1), bty="n")
163 7cf747b0 Jim Regetz
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")
164 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
165
166
plot(lats300, corByLat(d.bg, d.aster), type="l", xlab="Latitude",
167
    ylab="Correlation", ylim=ylim)
168
lines(lats150, corByLat(crop(d.bg, extent(d.srtm)), d.srtm), col="blue")
169
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
170
    lty=c(1, 1), bty="n")
171
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend")
172
abline(v=60, col="red", lty=2)
173
174
# ...with respect to CDEM
175
plot(lats300, corByLat(d.uncor, d.can), type="l", xlab="Latitude",
176
    ylab="Correlation", ylim=ylim)
177 7cf747b0 Jim Regetz
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fuse")
178 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
179
mtext(expression(paste(
180
    "Elevation correlations with respect to Canada DEM (",
181
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
182
183 7cf747b0 Jim Regetz
plot(lats300, corByLat(d.enblend, d.can), type="l", xlab="Latitude",
184 b02e39a0 Jim Regetz
    ylab="Correlation", ylim=ylim)
185 7cf747b0 Jim Regetz
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")
186 b02e39a0 Jim Regetz
abline(v=60, col="red", lty=2)
187
188
plot(lats300, corByLat(d.bg, d.can), type="l", xlab="Latitude",
189
    ylab="Correlation", ylim=ylim)
190
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend")
191
abline(v=60, col="red", lty=2)
192
193
# close pdf device driver
194
dev.off()
195
196 6626e7db Jim Regetz
#
197
# plot pattern of ASTER-SRTM deltas as a function of ASTER elevation
198
#
199
200
plotMeanDeltaByElev <- function(delta.vals, elev, ...) {
201
    mean.by.elev <- tapply(delta.vals, elev, mean)
202
    sd.by.elev <- tapply(delta.vals, elev, sd)
203
    n.by.elev <- tapply(delta.vals, elev, length)
204
    se.by.elev <- sd.by.elev/sqrt(n.by.elev)
205
    na.se.points <- mean.by.elev[is.na(se.by.elev)]
206
    se.by.elev[is.na(se.by.elev)] <- 0
207
    elev <- as.numeric(names(mean.by.elev))
208
    plot(elev, mean.by.elev, pch=16,
209
        xlim=c(0, max(elev)), ylim=c(min(mean.by.elev -
210
        se.by.elev), max(mean.by.elev + se.by.elev)), type="n", ...)
211
    segments(elev, mean.by.elev-se.by.elev,
212
        as.numeric(names(mean.by.elev)), mean.by.elev+se.by.elev,
213
        col="grey")
214
    points(elev, mean.by.elev, pch=".")
215
    points(as.numeric(names(na.se.points)), na.se.points, pch=4,
216
        col="red", cex=0.5)
217
}
218
219
220
d.aster.crop.vals <- values(crop(d.aster, extent(d.srtm)))
221
d.srtm.vals <- values(d.srtm)
222
delta.vals <- d.aster.crop.vals - d.srtm.vals
223
plotMeanDeltaByElev(delta.vals, d.aster.crop.vals,
224
    xlab="ASTER elevation (m)", ylab="ASTER-SRTM difference (m)")
225
226
plotDeltaBins <- function(delta.vals, elev, bin.min, bin.width, bin.max=1500,
227
    outline=FALSE, ...) {
228
    breaks <- seq(bin.min, bin.max, by=bin.width)
229
    midpts <- c(
230
        paste("<", bin.min, sep=""),
231
        head(breaks, -1) + bin.width/2,
232
        paste(">", bin.max, sep=""))
233
    elev <- cut(elev, breaks=c(0, breaks, Inf), labels=midpts)
234
    bp <- boxplot(delta.vals ~ elev, outline=outline, col="lightgray",
235
        frame=FALSE, ...)
236
    text(1:length(bp$n), bp$stats[5,], labels=round(bp$n/1000),
237
        pos=3, cex=0.5, offset=0.2, font=3, col="gray")
238
    #axis(3, at=seq_along(bp$n), labels=paste(round(bp$n/1000), "k", sep=""),
239
    #    cex.axis=0.7, tick=FALSE, font=3, line=-1)
240
    #mtext("n =", side=3, adj=0, font=3, cex=0.7)
241
    abline(h=median(delta.vals), col="red", lty=2)
242
    invisible(bp)
243
}
244
245
d.aster.crop.vals <- values(crop(d.aster, extent(d.srtm)))
246
d.srtm.vals <- values(d.srtm)
247
delta.vals <- d.aster.crop.vals - d.srtm.vals
248
#  d.aster.crop.vals <- d.aster.crop.vals[d.srtm.vals>0]
249
#  d.srtm.vals <- d.srtm.vals[d.srtm.vals>0]
250
251
png("aster-srtm-bins.png", height=5, width=8, units="in", res=300)
252 085ecdf3 Jim Regetz
plotDeltaBins(delta.vals, d.srtm.vals, 150, 50, 1500, las=2,
253
    cex.axis=0.8, xlab="Midpoints of SRTM elevation bins (m)",
254
    ylab="ASTER - SRTM difference (m)")
255 6626e7db Jim Regetz
dev.off()
256
257
# plot scatter of aster vs srtm
258
png("aster-srtm-scatter.png", height=5, width=8, units="in", res=300)
259
plot(jitter(d.srtm.vals), jitter(d.aster.crop.vals), pch=".",
260
    xlab="SRTM elevation (m)", ylab="ASTER elevation (m)", cex=0.5)
261
abline(median(delta.vals), 1, col="red", cex=0.5)
262
abline(0, 1, col="blue", lty=2, cex=0.5)
263 96629c37 Jim Regetz
# add inset histogram of differences
264
opar <- par(fig=c(0.55, 0.95, 0.1, 0.6), new=TRUE)
265
h <- hist(delta.vals[abs(delta.vals)<60], breaks=48, xlab=NA, main=NULL,
266
    col=grey(0.8), border=grey(0.3), yaxt="n", ylab=NA, cex.axis=0.5,
267
    cex.lab=0.5, tcl=-0.25, mgp=c(3,0,0))
268
text(10, 0.4*max(h$counts), labels=paste("Entire range:\n(",
269
    min(delta.vals), ", ", max(delta.vals), ")", sep=""), cex=0.6,
270
    adj=c(0,0))
271 8711607d Jim Regetz
mtext("ASTER - SRTM (m)", side=1, cex=0.5, line=0.6)
272 6626e7db Jim Regetz
dev.off()
273