Revision 5326e348
Added by Jim Regetz about 13 years ago
- ID 5326e348c68821d57764224c3dc24e809359d47c
| terrain/flow/flow-assessment.R | ||
|---|---|---|
| 1 |
# R code to plot latitudinal profiles of mean flow direction, along with
|
|
| 2 |
# both RMSE and correlation coefficients comparing fused layers with
|
|
| 3 |
# both the raw ASTER and with the Canada DEM |
|
| 1 |
# R code to plot latitudinal profiles of (circular) mean flow direction,
|
|
| 2 |
# along with both RMSE and (circular) correlation coefficients comparing
|
|
| 3 |
# fused layers with both the raw ASTER and with the Canada DEM
|
|
| 4 | 4 |
# |
| 5 | 5 |
# Jim Regetz |
| 6 | 6 |
# NCEAS |
| 7 |
# Created on 08-Jun-2011 |
|
| 8 | 7 |
|
| 9 | 8 |
library(raster) |
| 9 |
library(circular) |
|
| 10 | 10 |
|
| 11 | 11 |
datadir <- "/home/regetz/media/temp/terrain/flow" |
| 12 | 12 |
|
| 13 |
# create function to recode values into degrees |
|
| 13 |
# create function to recode terraflow SFD values into degrees, with |
|
| 14 |
# 0=North and proceeding clockwise (this matches gdaldem's default |
|
| 15 |
# azimuth output for aspect calculation) |
|
| 14 | 16 |
recode <- function(r) {
|
| 15 | 17 |
v <- values(r) |
| 16 | 18 |
v[v==0] <- NA |
| 17 |
v[v==1] <- 0 |
|
| 18 |
v[v==2] <- 45 |
|
| 19 |
v[v==3] <- 90 |
|
| 20 |
v[v==4] <- 90 |
|
| 21 |
v[v==8] <- 135 |
|
| 22 |
v[v==16] <- 180 |
|
| 23 |
v[v==32] <- 225 |
|
| 24 |
v[v==64] <- 270 |
|
| 25 |
v[v==128] <- 315 |
|
| 19 |
v[v==1] <- 90 ## east |
|
| 20 |
v[v==2] <- 135 |
|
| 21 |
v[v==4] <- 180 ## south |
|
| 22 |
v[v==8] <- 225 |
|
| 23 |
v[v==16] <- 270 ## west |
|
| 24 |
v[v==32] <- 315 |
|
| 25 |
v[v==64] <- 0 ## north |
|
| 26 |
v[v==128] <- 45 |
|
| 26 | 27 |
r[] <- v |
| 27 | 28 |
return(r) |
| 28 | 29 |
} |
| ... | ... | |
| 31 | 32 |
sfd.aster <- recode(raster(file.path(datadir, "aster_300straddle_sfd.tif"))) |
| 32 | 33 |
sfd.srtm <- recode(raster(file.path(datadir, "srtm_150below_sfd.tif"))) |
| 33 | 34 |
sfd.uncor <- recode(raster(file.path(datadir, "fused_300straddle_sfd.tif"))) |
| 34 |
#sfd.eramp <- recode(raster(file.path(datadir, "fused_300straddle_rampexp_sfd.tif")))
|
|
| 35 |
sfd.enblend <- recode(raster(file.path(datadir, "fused_300straddle_enblend_sfd.tif")))
|
|
| 35 | 36 |
sfd.bg <- recode(raster(file.path(datadir, "fused_300straddle_blendgau_sfd.tif"))) |
| 36 | 37 |
sfd.can <- recode(raster(file.path(datadir, "cdem_300straddle_sfd.tif"))) |
| 37 | 38 |
|
| ... | ... | |
| 46 | 47 |
# plot latitudinal profiles of mean flow direction |
| 47 | 48 |
# |
| 48 | 49 |
|
| 49 |
par(mfrow=c(2,2), omi=c(1,1,1,1)) |
|
| 50 |
# simple helper function to calculate row-wise means using circular |
|
| 51 |
# mean, patterned after circ.mean in the CircStats package |
|
| 52 |
rowMeansC <- function(r1, na.rm=TRUE) {
|
|
| 53 |
m1 <- as.matrix(r1) |
|
| 54 |
m1[] <- (m1 * pi)/180 |
|
| 55 |
sinr <- rowSums(sin(m1), na.rm=na.rm) |
|
| 56 |
cosr <- rowSums(cos(m1), na.rm=na.rm) |
|
| 57 |
cmeans <- atan2(sinr, cosr) |
|
| 58 |
(cmeans * 180)/pi |
|
| 59 |
} |
|
| 50 | 60 |
|
| 51 |
ylim <- c(80, 280) |
|
| 61 |
par(mfrow=c(2,2), omi=c(1,1,1,1)) |
|
| 62 |
ylim <- c(-180, 180) |
|
| 52 | 63 |
|
| 53 |
plot(lats300, rowMeans(as.matrix(sfd.uncor), na.rm=TRUE), type="l",
|
|
| 64 |
plot(lats300, rowMeansC(sfd.can), type="l", yaxt="n",
|
|
| 54 | 65 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
| 55 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="uncorrected") |
|
| 66 |
axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))
|
|
| 67 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Original DEMs") |
|
| 68 |
lines(lats300, rowMeansC(sfd.aster), col="blue") |
|
| 69 |
lines(lats150, rowMeansC(sfd.srtm), col="red") |
|
| 70 |
legend("bottomright", legend=c("ASTER", "SRTM", "CDED"), col=c("blue",
|
|
| 71 |
"red", "black"), lty=c(1, 1), bty="n") |
|
| 56 | 72 |
abline(v=60, col="red", lty=2) |
| 57 | 73 |
mtext(expression(paste("Latitudinal profiles of mean flow direction (",
|
| 58 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
|
|
| 74 |
125*degree, "W to ", 100*degree, "W)")), adj=0, line=2, font=2)
|
|
| 59 | 75 |
|
| 60 |
plot(lats300, rowMeans(as.matrix(sfd.can), na.rm=TRUE), type="l", |
|
| 76 |
#plot(lats300, rowMeans(as.matrix(sfd.uncor), na.rm=TRUE), type="l", |
|
| 77 |
# xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
|
| 78 |
#text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="uncorrected") |
|
| 79 |
#abline(v=60, col="red", lty=2) |
|
| 80 |
#mtext(expression(paste("Latitudinal profiles of mean flow direction (",
|
|
| 81 |
# 125*degree, "W to ", 100*degree, "W)")), adj=0, line=2, font=2) |
|
| 82 |
|
|
| 83 |
plot(lats300, rowMeansC(sfd.uncor), type="l", yaxt="n", |
|
| 61 | 84 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
| 62 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Canada DEM") |
|
| 85 |
axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))
|
|
| 86 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="simple fused") |
|
| 63 | 87 |
abline(v=60, col="red", lty=2) |
| 64 | 88 |
|
| 65 |
#plot(lats300, rowMeans(as.matrix(sfd.eramp), na.rm=TRUE), type="l", |
|
| 66 |
plot(lats300, rowMeans(as.matrix(sfd.can), na.rm=TRUE), type="n", |
|
| 89 |
plot(lats300, rowMeansC(sfd.enblend), type="l", yaxt="n", |
|
| 67 | 90 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
| 68 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="exponential ramp")
|
|
| 69 |
text(mean(lats300), mean(ylim), pos=1, font=3, labels="(skipped)")
|
|
| 70 |
#abline(v=60, col="red", lty=2)
|
|
| 91 |
axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))
|
|
| 92 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="multires spline")
|
|
| 93 |
abline(v=60, col="red", lty=2) |
|
| 71 | 94 |
|
| 72 |
plot(lats300, rowMeans(as.matrix(sfd.bg), na.rm=TRUE), type="l",
|
|
| 95 |
plot(lats300, rowMeansC(sfd.bg), type="l", yaxt="n",
|
|
| 73 | 96 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
| 97 |
axis(2, at=c(-180, -90, 0, 90, 180), labels=c("S", "W", "N", "E", "S"))
|
|
| 74 | 98 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend") |
| 75 | 99 |
abline(v=60, col="red", lty=2) |
| 76 | 100 |
|
| 77 | 101 |
|
| 102 |
|
|
| 78 | 103 |
# |
| 79 | 104 |
# plot latitudinal profiles of RMSE |
| 80 | 105 |
# |
| ... | ... | |
| 103 | 128 |
abline(v=60, col="red", lty=2) |
| 104 | 129 |
mtext(expression(paste( |
| 105 | 130 |
"Flowdir discrepancies with respect to separate ASTER/SRTM components (",
|
| 106 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
|
|
| 131 |
125*degree, "W to ", 100*degree, "W)")), adj=0, line=2, font=2)
|
|
| 107 | 132 |
|
| 108 |
#plot(lats300, rmse(sfd.eramp, sfd.aster), type="l", xlab="Latitude", |
|
| 109 |
plot(lats300, rep(0, 300), type="n", xlab="Latitude", |
|
| 133 |
plot(lats300, rmse(sfd.enblend, sfd.aster), type="l", xlab="Latitude", |
|
| 110 | 134 |
ylab="RMSE", ylim=ylim) |
| 111 |
#lines(lats150, rmse(crop(sfd.eramp, extent(sfd.srtm)), sfd.srtm), col="blue")
|
|
| 112 |
#legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 113 |
# lty=c(1, 1), bty="n")
|
|
| 135 |
lines(lats150, rmse(crop(sfd.enblend, extent(sfd.srtm)), sfd.srtm), col="blue")
|
|
| 136 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 137 |
lty=c(1, 1), bty="n") |
|
| 114 | 138 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="exponential ramp") |
| 115 | 139 |
text(mean(lats300), mean(ylim), font=3, labels="(skipped)") |
| 116 |
#abline(v=60, col="red", lty=2)
|
|
| 140 |
abline(v=60, col="red", lty=2) |
|
| 117 | 141 |
|
| 118 | 142 |
plot(lats300, rmse(sfd.bg, sfd.aster), type="l", xlab="Latitude", |
| 119 | 143 |
ylab="RMSE", ylim=ylim) |
| ... | ... | |
| 130 | 154 |
abline(v=60, col="red", lty=2) |
| 131 | 155 |
mtext(expression(paste( |
| 132 | 156 |
"Flowdir discrepancies with respect to Canada DEM (",
|
| 133 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2)
|
|
| 157 |
125*degree, "W to ", 100*degree, "W)")), adj=0, line=2, font=2)
|
|
| 134 | 158 |
|
| 135 |
#plot(lats300, rmse(sfd.eramp, sfd.can), type="l", xlab="Latitude", |
|
| 136 |
plot(lats300, rep(0, 300), type="n", xlab="Latitude", |
|
| 159 |
plot(lats300, rmse(sfd.enblend, sfd.can), type="l", xlab="Latitude", |
|
| 137 | 160 |
ylab="RMSE", ylim=ylim) |
| 138 | 161 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="exponential ramp") |
| 139 | 162 |
text(mean(lats300), mean(ylim), font=3, labels="(skipped)") |
| 140 |
#abline(v=60, col="red", lty=2)
|
|
| 163 |
abline(v=60, col="red", lty=2) |
|
| 141 | 164 |
|
| 142 | 165 |
plot(lats300, rmse(sfd.bg, sfd.can), type="l", xlab="Latitude", |
| 143 | 166 |
ylab="RMSE", ylim=ylim) |
| 144 | 167 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend") |
| 145 | 168 |
abline(v=60, col="red", lty=2) |
| 146 | 169 |
|
| 170 |
|
|
| 171 |
# |
|
| 172 |
# plot latitudinal profiles of correlation coefficients |
|
| 173 |
# |
|
| 174 |
|
|
| 175 |
# simple helper function to calculate row-wise *circular* correlation |
|
| 176 |
# coefficients |
|
| 177 |
corByLat <- function(r1, r2, rows) {
|
|
| 178 |
if (missing(rows)) {
|
|
| 179 |
rows <- 1:nrow(r1) |
|
| 180 |
} |
|
| 181 |
m1 <- circular(as.matrix(r1), units="degrees", rotation="clock") |
|
| 182 |
m2 <- circular(as.matrix(r2), units="degrees", rotation="clock") |
|
| 183 |
sapply(rows, function(row) {
|
|
| 184 |
p <- cor.circular(m1[row,], m2[row,]) |
|
| 185 |
if (is.null(p)) NA else p |
|
| 186 |
}) |
|
| 187 |
} |
|
| 188 |
|
|
| 189 |
par(mfrow=c(2,3), omi=c(1,1,1,1)) |
|
| 190 |
|
|
| 191 |
ylim <- c(-1, 1) |
|
| 192 |
|
|
| 193 |
# ...with respect to ASTER |
|
| 194 |
plot(lats300, corByLat(sfd.uncor, sfd.aster), type="l", xlab="Latitude", |
|
| 195 |
ylab="Circular correlation", ylim=ylim) |
|
| 196 |
lines(lats150, corByLat(crop(sfd.uncor, extent(sfd.srtm)), sfd.srtm), col="blue") |
|
| 197 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 198 |
lty=c(1, 1), bty="n") |
|
| 199 |
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused") |
|
| 200 |
abline(v=60, col="red", lty=2) |
|
| 201 |
mtext(expression(paste( |
|
| 202 |
"Flow direction correlations with respect to separate ASTER/SRTM components (",
|
|
| 203 |
125*degree, "W to ", 100*degree, "W)")), adj=0, line=2, font=2) |
|
| 204 |
|
|
| 205 |
plot(lats300, corByLat(sfd.enblend, sfd.aster), type="l", xlab="Latitude", |
|
| 206 |
ylab="Circular correlation", ylim=ylim) |
|
| 207 |
lines(lats150, corByLat(crop(sfd.enblend, extent(sfd.srtm)), sfd.srtm), col="blue") |
|
| 208 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 209 |
lty=c(1, 1), bty="n") |
|
| 210 |
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline") |
|
| 211 |
abline(v=60, col="red", lty=2) |
|
| 212 |
|
|
| 213 |
plot(lats300, corByLat(sfd.bg, sfd.aster), type="l", xlab="Latitude", |
|
| 214 |
ylab="Circular correlation", ylim=ylim) |
|
| 215 |
lines(lats150, corByLat(crop(sfd.bg, extent(sfd.srtm)), sfd.srtm), col="blue") |
|
| 216 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 217 |
lty=c(1, 1), bty="n") |
|
| 218 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
|
| 219 |
abline(v=60, col="red", lty=2) |
|
| 220 |
|
|
| 221 |
# ...with respect to CDEM |
|
| 222 |
plot(lats300, corByLat(sfd.uncor, sfd.can), type="l", xlab="Latitude", |
|
| 223 |
ylab="Circular correlation", ylim=ylim) |
|
| 224 |
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused") |
|
| 225 |
abline(v=60, col="red", lty=2) |
|
| 226 |
mtext(expression(paste( |
|
| 227 |
"Flow direction correlations with respect to Canada DEM (",
|
|
| 228 |
125*degree, "W to ", 100*degree, "W)")), adj=0, line=2, font=2) |
|
| 229 |
|
|
| 230 |
plot(lats300, corByLat(sfd.enblend, sfd.can), type="l", xlab="Latitude", |
|
| 231 |
ylab="Circular correlation", ylim=ylim) |
|
| 232 |
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline") |
|
| 233 |
abline(v=60, col="red", lty=2) |
|
| 234 |
|
|
| 235 |
plot(lats300, corByLat(sfd.bg, sfd.can), type="l", xlab="Latitude", |
|
| 236 |
ylab="Circular correlation", ylim=ylim) |
|
| 237 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
|
| 238 |
abline(v=60, col="red", lty=2) |
|
| 239 |
|
|
| 147 | 240 |
# close pdf device driver |
| 148 | 241 |
dev.off() |
Also available in: Unified diff
revamped aspect/flowdir boundary assessments to use circular statistics