Project

General

Profile

Download (7.03 KB) Statistics
| Branch: | Revision:
1
# R code to plot latitudinal profiles of mean slope, 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/slope"
12

    
13
# load slope rasters
14
s.aster <- raster(file.path(datadir, "aster_300straddle_s.tif"))
15
s.srtm <- raster(file.path(datadir, "srtm_150below_s.tif"))
16
s.uncor <- raster(file.path(datadir, "fused_300straddle_s.tif"))
17
s.enblend <- raster(file.path(datadir, "fused_300straddle_enblend_s.tif"))
18
s.bg <- raster(file.path(datadir, "fused_300straddle_blendgau_s.tif"))
19
s.can <- raster(file.path(datadir, "cdem_300straddle_s.tif"))
20

    
21
# extract raster latitudes for later
22
lats300 <- yFromRow(s.aster, 1:nrow(s.aster))
23
lats150 <- yFromRow(s.srtm, 1:nrow(s.srtm))
24

    
25
# initialize output pdf device driver
26
pdf("slope-assessment.pdf", height=8, width=11.5)
27

    
28

    
29
#
30
# plot latitudinal profiles of mean slope
31
#
32

    
33
par(mfrow=c(2,2), omi=c(1,1,1,1))
34

    
35
ylim <- c(1, 6)
36

    
37
plot(lats300, rowMeans(as.matrix(s.can), na.rm=TRUE), type="l",
38
    xlab="Latitude", ylab="Mean slope", ylim=ylim)
39
text(min(lats300), max(ylim)-0.5, pos=4, font=3, labels="Original DEMs")
40
lines(lats300, rowMeans(as.matrix(s.aster), na.rm=TRUE), col="blue")
41
lines(lats150, rowMeans(as.matrix(s.srtm), na.rm=TRUE), col="red")
42
legend("topright", legend=c("ASTER", "SRTM", "CDED"), col=c("blue",
43
    "red", "black"), lty=c(1, 1), bty="n")
44
abline(v=60, col="red", lty=2)
45
mtext(expression(paste("Latitudinal profiles of mean slope (",
46
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
47

    
48
plot(lats300, rowMeans(as.matrix(s.uncor), na.rm=TRUE), type="l",
49
    xlab="Latitude", ylab="Mean slope", ylim=ylim)
50
text(min(lats300), max(ylim)-0.5, pos=4, font=3, labels="simple fused")
51
abline(v=60, col="red", lty=2)
52

    
53
plot(lats300, rowMeans(as.matrix(s.enblend), na.rm=TRUE), type="l",
54
    xlab="Latitude", ylab="Mean slope", ylim=ylim)
55
text(min(lats300), max(ylim)-0.5, pos=4, font=3, labels="multires spline")
56
abline(v=60, col="red", lty=2)
57

    
58
plot(lats300, rowMeans(as.matrix(s.bg), na.rm=TRUE), type="l",
59
    xlab="Latitude", ylab="Mean slope", ylim=ylim)
60
text(min(lats300), max(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
# ...with respect to ASTER
78
plot(lats300, rmse(s.uncor, s.aster), type="l", xlab="Latitude",
79
    ylab="RMSE", ylim=c(0, 5))
80
lines(lats150, rmse(crop(s.uncor, extent(s.srtm)), s.srtm), col="blue")
81
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
82
    lty=c(1, 1), bty="n")
83
text(min(lats300), 4.5, pos=4, font=3, labels="simple fused")
84
abline(v=60, col="red", lty=2)
85
mtext(expression(paste(
86
    "Slope discrepancies with respect to separate ASTER/SRTM components (",
87
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
88

    
89
plot(lats300, rmse(s.enblend, s.aster), type="l", xlab="Latitude",
90
    ylab="RMSE", ylim=c(0, 5))
91
lines(lats150, rmse(crop(s.enblend, extent(s.srtm)), s.srtm), col="blue")
92
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
93
    lty=c(1, 1), bty="n")
94
text(min(lats300), 4.5, pos=4, font=3, labels="multires spline")
95
abline(v=60, col="red", lty=2)
96

    
97
plot(lats300, rmse(s.bg, s.aster), type="l", xlab="Latitude",
98
    ylab="RMSE", ylim=c(0, 5))
99
lines(lats150, rmse(crop(s.bg, extent(s.srtm)), s.srtm), col="blue")
100
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
101
    lty=c(1, 1), bty="n")
102
text(min(lats300), 4.5, pos=4, font=3, labels="gaussian blend")
103
abline(v=60, col="red", lty=2)
104

    
105
# ...with respect to CDEM
106
plot(lats300, rmse(s.uncor, s.can), type="l", xlab="Latitude",
107
    ylab="RMSE", ylim=c(0, 5))
108
text(min(lats300), 4.5, pos=4, font=3, labels="simple fused")
109
abline(v=60, col="red", lty=2)
110
mtext(expression(paste(
111
    "Slope discrepancies with respect to Canada DEM (",
112
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
113

    
114
plot(lats300, rmse(s.enblend, s.can), type="l", xlab="Latitude",
115
    ylab="RMSE", ylim=c(0, 5))
116
text(min(lats300), 4.5, pos=4, font=3, labels="multires spline")
117
abline(v=60, col="red", lty=2)
118

    
119
plot(lats300, rmse(s.bg, s.can), type="l", xlab="Latitude",
120
    ylab="RMSE", ylim=c(0, 5))
121
text(min(lats300), 4.5, pos=4, font=3, labels="gaussian blend")
122
abline(v=60, col="red", lty=2)
123

    
124

    
125
#
126
# plot latitudinal profiles of correlation coefficients
127
#
128

    
129
# simple helper function to calculate row-wise correlation coefficients
130
corByLat <- function(r1, r2, rows) {
131
    if (missing(rows)) {
132
        rows <- 1:nrow(r1)
133
    }
134
    m1 <- as.matrix(r1)
135
    m2 <- as.matrix(r2)
136
    sapply(rows, function(row) cor(m1[row,], m2[row,],
137
        use="pairwise.complete.obs"))
138
}
139

    
140
par(mfrow=c(2,3), omi=c(1,1,1,1))
141

    
142
ylim <- c(0.65, 1)
143

    
144
# ...with respect to ASTER
145
plot(lats300, corByLat(s.uncor, s.aster), type="l", xlab="Latitude",
146
    ylab="Correlation", ylim=ylim)
147
lines(lats150, corByLat(crop(s.uncor, extent(s.srtm)), s.srtm), col="blue")
148
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
149
    lty=c(1, 1), bty="n")
150
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused")
151
abline(v=60, col="red", lty=2)
152
mtext(expression(paste(
153
    "Slope correlations with respect to separate ASTER/SRTM components (",
154
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
155

    
156
plot(lats300, corByLat(s.enblend, s.aster), type="l", xlab="Latitude",
157
    ylab="Correlation", ylim=ylim)
158
lines(lats150, corByLat(crop(s.enblend, extent(s.srtm)), s.srtm), col="blue")
159
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
160
    lty=c(1, 1), bty="n")
161
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")
162
abline(v=60, col="red", lty=2)
163

    
164
plot(lats300, corByLat(s.bg, s.aster), type="l", xlab="Latitude",
165
    ylab="Correlation", ylim=ylim)
166
lines(lats150, corByLat(crop(s.bg, extent(s.srtm)), s.srtm), col="blue")
167
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
168
    lty=c(1, 1), bty="n")
169
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend")
170
abline(v=60, col="red", lty=2)
171

    
172
# ...with respect to CDEM
173
plot(lats300, corByLat(s.uncor, s.can), type="l", xlab="Latitude",
174
    ylab="Correlation", ylim=ylim)
175
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused")
176
abline(v=60, col="red", lty=2)
177
mtext(expression(paste(
178
    "Slope correlations with respect to Canada DEM (",
179
    136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
180

    
181
plot(lats300, corByLat(s.enblend, s.can), type="l", xlab="Latitude",
182
    ylab="Correlation", ylim=ylim)
183
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")
184
abline(v=60, col="red", lty=2)
185

    
186
plot(lats300, corByLat(s.bg, s.can), type="l", xlab="Latitude",
187
    ylab="Correlation", ylim=ylim)
188
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend")
189
abline(v=60, col="red", lty=2)
190

    
191
# close pdf device driver
192
dev.off()
    (1-1/1)