Revision b02e39a0
Added by Jim Regetz about 13 years ago
- ID b02e39a0b4de755bb216f9272065a6d838530cc1
| terrain/aspect/aspect-assessment.R | ||
|---|---|---|
| 61 | 61 |
# plot latitudinal profiles of RMSE |
| 62 | 62 |
# |
| 63 | 63 |
|
| 64 |
# simple helper function to calculate row-wise RMSEs |
|
| 65 |
rmse <- function(r1, r2, na.rm=TRUE) {
|
|
| 66 |
sqrt(rowMeans(as.matrix((r1 - r2)^2), na.rm=na.rm)) |
|
| 64 |
# simple helper function to calculate row-wise RMSEs, accounting for the |
|
| 65 |
# fact that aspect values are circular (0-360), so the difference |
|
| 66 |
# between e.g. 5 and 355 should only be 10 |
|
| 67 |
rmse <- function(r1, r2, na.rm=TRUE, use) {
|
|
| 68 |
diffs <- abs(as.matrix(r1) - as.matrix(r2)) |
|
| 69 |
if (!missing(use)) diffs[!use] <- NA |
|
| 70 |
diffs[] <- ifelse(diffs>180, 360-diffs, diffs) |
|
| 71 |
sqrt(rowMeans(diffs^2, na.rm=na.rm)) |
|
| 67 | 72 |
} |
| 68 | 73 |
|
| 69 | 74 |
par(mfrow=c(2,3), omi=c(1,1,1,1)) |
| 70 | 75 |
|
| 71 |
ylim <- c(0, 170)
|
|
| 76 |
ylim <- c(0, 100)
|
|
| 72 | 77 |
|
| 73 | 78 |
# ...with respect to ASTER |
| 74 | 79 |
plot(lats300, rmse(s.uncor, s.aster), type="l", xlab="Latitude", |
| ... | ... | |
| 117 | 122 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend") |
| 118 | 123 |
abline(v=60, col="red", lty=2) |
| 119 | 124 |
|
| 125 |
# close pdf device driver |
|
| 126 |
dev.off() |
|
| 127 |
|
|
| 128 |
stop("not doing correlations")
|
|
| 120 | 129 |
|
| 121 | 130 |
# |
| 122 | 131 |
# plot latitudinal profiles of correlation coefficients |
| ... | ... | |
| 183 | 192 |
ylab="Correlation", ylim=ylim) |
| 184 | 193 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
| 185 | 194 |
abline(v=60, col="red", lty=2) |
| 186 |
|
|
| 187 |
# close pdf device driver |
|
| 188 |
dev.off() |
|
| terrain/dem/dem-assessment.R | ||
|---|---|---|
| 1 |
# 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 |
d.eramp <- raster(file.path(datadir, "fused_300straddle_rampexp.tif")) |
|
| 18 |
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 |
# initialize output pdf device driver |
|
| 26 |
pdf("elevation-assessment.pdf", height=8, width=11.5)
|
|
| 27 |
|
|
| 28 |
|
|
| 29 |
# |
|
| 30 |
# plot latitudinal profiles of mean elevation |
|
| 31 |
# |
|
| 32 |
|
|
| 33 |
par(mfrow=c(2,2), omi=c(1,1,1,1)) |
|
| 34 |
|
|
| 35 |
ylim <- c(550, 575) |
|
| 36 |
|
|
| 37 |
plot(lats300, rowMeans(as.matrix(d.uncor), na.rm=TRUE), type="l", |
|
| 38 |
xlab="Latitude", ylab="Mean elevation", ylim=ylim) |
|
| 39 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="uncorrected") |
|
| 40 |
abline(v=60, col="red", lty=2) |
|
| 41 |
mtext(expression(paste("Latitudinal profiles of mean elevation (",
|
|
| 42 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 43 |
|
|
| 44 |
plot(lats300, rowMeans(as.matrix(d.can), na.rm=TRUE), type="l", |
|
| 45 |
xlab="Latitude", ylab="Mean elevation", ylim=ylim) |
|
| 46 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Canada DEM") |
|
| 47 |
abline(v=60, col="red", lty=2) |
|
| 48 |
|
|
| 49 |
plot(lats300, rowMeans(as.matrix(d.eramp), na.rm=TRUE), type="l", |
|
| 50 |
xlab="Latitude", ylab="Mean elevation", ylim=ylim) |
|
| 51 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="exponential ramp") |
|
| 52 |
abline(v=60, col="red", lty=2) |
|
| 53 |
|
|
| 54 |
plot(lats300, rowMeans(as.matrix(d.bg), na.rm=TRUE), type="l", |
|
| 55 |
xlab="Latitude", ylab="Mean elevation", ylim=ylim) |
|
| 56 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend") |
|
| 57 |
abline(v=60, col="red", lty=2) |
|
| 58 |
|
|
| 59 |
|
|
| 60 |
# |
|
| 61 |
# plot latitudinal profiles of RMSE |
|
| 62 |
# |
|
| 63 |
|
|
| 64 |
# simple helper function to calculate row-wise RMSEs |
|
| 65 |
rmse <- function(r1, r2, na.rm=TRUE, use) {
|
|
| 66 |
diffs <- abs(as.matrix(r1) - as.matrix(r2)) |
|
| 67 |
if (!missing(use)) diffs[!use] <- NA |
|
| 68 |
sqrt(rowMeans(diffs^2, na.rm=na.rm)) |
|
| 69 |
} |
|
| 70 |
|
|
| 71 |
par(mfrow=c(2,3), omi=c(1,1,1,1)) |
|
| 72 |
|
|
| 73 |
ylim <- c(0, 35) |
|
| 74 |
|
|
| 75 |
# ...with respect to ASTER |
|
| 76 |
plot(lats300, rmse(d.uncor, d.aster), type="l", xlab="Latitude", |
|
| 77 |
ylab="RMSE", ylim=ylim) |
|
| 78 |
lines(lats150, rmse(crop(d.uncor, extent(d.srtm)), d.srtm), col="blue") |
|
| 79 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 80 |
lty=c(1, 1), bty="n") |
|
| 81 |
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="uncorrected") |
|
| 82 |
abline(v=60, col="red", lty=2) |
|
| 83 |
mtext(expression(paste( |
|
| 84 |
"Elevation discrepancies with respect to separate ASTER/SRTM components (",
|
|
| 85 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 86 |
|
|
| 87 |
plot(lats300, rmse(d.eramp, d.aster), type="l", xlab="Latitude", |
|
| 88 |
ylab="RMSE", ylim=ylim) |
|
| 89 |
lines(lats150, rmse(crop(d.eramp, extent(d.srtm)), d.srtm), col="blue") |
|
| 90 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 91 |
lty=c(1, 1), bty="n") |
|
| 92 |
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="exponential ramp") |
|
| 93 |
abline(v=60, col="red", lty=2) |
|
| 94 |
|
|
| 95 |
plot(lats300, rmse(d.bg, d.aster), type="l", xlab="Latitude", |
|
| 96 |
ylab="RMSE", ylim=ylim) |
|
| 97 |
lines(lats150, rmse(crop(d.bg, extent(d.srtm)), d.srtm), col="blue") |
|
| 98 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 99 |
lty=c(1, 1), bty="n") |
|
| 100 |
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="gaussian blend") |
|
| 101 |
abline(v=60, col="red", lty=2) |
|
| 102 |
|
|
| 103 |
# ...with respect to CDEM |
|
| 104 |
plot(lats300, rmse(d.uncor, d.can), type="l", xlab="Latitude", |
|
| 105 |
ylab="RMSE", ylim=ylim) |
|
| 106 |
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="uncorrected") |
|
| 107 |
abline(v=60, col="red", lty=2) |
|
| 108 |
mtext(expression(paste( |
|
| 109 |
"Elevation discrepancies with respect to Canada DEM (",
|
|
| 110 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 111 |
|
|
| 112 |
plot(lats300, rmse(d.eramp, d.can), type="l", xlab="Latitude", |
|
| 113 |
ylab="RMSE", ylim=ylim) |
|
| 114 |
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="exponential ramp") |
|
| 115 |
abline(v=60, col="red", lty=2) |
|
| 116 |
|
|
| 117 |
plot(lats300, rmse(d.bg, d.can), type="l", xlab="Latitude", |
|
| 118 |
ylab="RMSE", ylim=ylim) |
|
| 119 |
text(min(lats300), max(ylim)-1, pos=4, font=3, labels="gaussian blend") |
|
| 120 |
abline(v=60, col="red", lty=2) |
|
| 121 |
|
|
| 122 |
|
|
| 123 |
# |
|
| 124 |
# plot latitudinal profiles of correlation coefficients |
|
| 125 |
# |
|
| 126 |
|
|
| 127 |
# simple helper function to calculate row-wise correlation coefficients |
|
| 128 |
corByLat <- function(r1, r2, rows) {
|
|
| 129 |
if (missing(rows)) {
|
|
| 130 |
rows <- 1:nrow(r1) |
|
| 131 |
} |
|
| 132 |
m1 <- as.matrix(r1) |
|
| 133 |
m2 <- as.matrix(r2) |
|
| 134 |
sapply(rows, function(row) cor(m1[row,], m2[row,], |
|
| 135 |
use="pairwise.complete.obs")) |
|
| 136 |
} |
|
| 137 |
|
|
| 138 |
par(mfrow=c(2,3), omi=c(1,1,1,1)) |
|
| 139 |
|
|
| 140 |
ylim <- c(0.99, 1) |
|
| 141 |
|
|
| 142 |
# ...with respect to ASTER |
|
| 143 |
plot(lats300, corByLat(d.uncor, d.aster), type="l", xlab="Latitude", |
|
| 144 |
ylab="Correlation", ylim=ylim) |
|
| 145 |
lines(lats150, corByLat(crop(d.uncor, extent(d.srtm)), d.srtm), col="blue") |
|
| 146 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 147 |
lty=c(1, 1), bty="n") |
|
| 148 |
text(min(lats300), min(ylim), pos=4, font=3, labels="uncorrected") |
|
| 149 |
abline(v=60, col="red", lty=2) |
|
| 150 |
mtext(expression(paste( |
|
| 151 |
"Elevation correlations with respect to separate ASTER/SRTM components (",
|
|
| 152 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 153 |
|
|
| 154 |
plot(lats300, corByLat(d.eramp, d.aster), type="l", xlab="Latitude", |
|
| 155 |
ylab="Correlation", ylim=ylim) |
|
| 156 |
lines(lats150, corByLat(crop(d.eramp, extent(d.srtm)), d.srtm), col="blue") |
|
| 157 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 158 |
lty=c(1, 1), bty="n") |
|
| 159 |
text(min(lats300), min(ylim), pos=4, font=3, labels="exponential ramp") |
|
| 160 |
abline(v=60, col="red", lty=2) |
|
| 161 |
|
|
| 162 |
plot(lats300, corByLat(d.bg, d.aster), type="l", xlab="Latitude", |
|
| 163 |
ylab="Correlation", ylim=ylim) |
|
| 164 |
lines(lats150, corByLat(crop(d.bg, extent(d.srtm)), d.srtm), col="blue") |
|
| 165 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 166 |
lty=c(1, 1), bty="n") |
|
| 167 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
|
| 168 |
abline(v=60, col="red", lty=2) |
|
| 169 |
|
|
| 170 |
# ...with respect to CDEM |
|
| 171 |
plot(lats300, corByLat(d.uncor, d.can), type="l", xlab="Latitude", |
|
| 172 |
ylab="Correlation", ylim=ylim) |
|
| 173 |
text(min(lats300), min(ylim), pos=4, font=3, labels="uncorrected") |
|
| 174 |
abline(v=60, col="red", lty=2) |
|
| 175 |
mtext(expression(paste( |
|
| 176 |
"Elevation correlations with respect to Canada DEM (",
|
|
| 177 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 178 |
|
|
| 179 |
plot(lats300, corByLat(d.eramp, d.can), type="l", xlab="Latitude", |
|
| 180 |
ylab="Correlation", ylim=ylim) |
|
| 181 |
text(min(lats300), min(ylim), pos=4, font=3, labels="exponential ramp") |
|
| 182 |
abline(v=60, col="red", lty=2) |
|
| 183 |
|
|
| 184 |
plot(lats300, corByLat(d.bg, d.can), type="l", xlab="Latitude", |
|
| 185 |
ylab="Correlation", ylim=ylim) |
|
| 186 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
|
| 187 |
abline(v=60, col="red", lty=2) |
|
| 188 |
|
|
| 189 |
# close pdf device driver |
|
| 190 |
dev.off() |
|
| 191 |
|
|
| 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 |
|
| 4 |
# |
|
| 5 |
# Jim Regetz |
|
| 6 |
# NCEAS |
|
| 7 |
# Created on 08-Jun-2011 |
|
| 8 |
|
|
| 9 |
library(raster) |
|
| 10 |
|
|
| 11 |
datadir <- "/home/regetz/media/temp/terrain/flow" |
|
| 12 |
|
|
| 13 |
# create function to recode values into degrees |
|
| 14 |
recode <- function(r) {
|
|
| 15 |
v <- values(r) |
|
| 16 |
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 |
|
| 26 |
r[] <- v |
|
| 27 |
return(r) |
|
| 28 |
} |
|
| 29 |
|
|
| 30 |
# load flow direction rasters, recoding on the fly |
|
| 31 |
sfd.aster <- recode(raster(file.path(datadir, "aster_300straddle_sfd.tif"))) |
|
| 32 |
sfd.srtm <- recode(raster(file.path(datadir, "srtm_150below_sfd.tif"))) |
|
| 33 |
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.bg <- recode(raster(file.path(datadir, "fused_300straddle_blendgau_sfd.tif"))) |
|
| 36 |
sfd.can <- recode(raster(file.path(datadir, "cdem_300straddle_sfd.tif"))) |
|
| 37 |
|
|
| 38 |
# extract raster latitudes for later |
|
| 39 |
lats300 <- yFromRow(sfd.aster, 1:nrow(sfd.aster)) |
|
| 40 |
lats150 <- yFromRow(sfd.srtm, 1:nrow(sfd.srtm)) |
|
| 41 |
|
|
| 42 |
# initialize output pdf device driver |
|
| 43 |
pdf("flowdir-assessment.pdf", height=8, width=11.5)
|
|
| 44 |
|
|
| 45 |
# |
|
| 46 |
# plot latitudinal profiles of mean flow direction |
|
| 47 |
# |
|
| 48 |
|
|
| 49 |
par(mfrow=c(2,2), omi=c(1,1,1,1)) |
|
| 50 |
|
|
| 51 |
ylim <- c(80, 280) |
|
| 52 |
|
|
| 53 |
plot(lats300, rowMeans(as.matrix(sfd.uncor), na.rm=TRUE), type="l", |
|
| 54 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
|
| 55 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="uncorrected") |
|
| 56 |
abline(v=60, col="red", lty=2) |
|
| 57 |
mtext(expression(paste("Latitudinal profiles of mean flow direction (",
|
|
| 58 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 59 |
|
|
| 60 |
plot(lats300, rowMeans(as.matrix(sfd.can), na.rm=TRUE), type="l", |
|
| 61 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
|
| 62 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Canada DEM") |
|
| 63 |
abline(v=60, col="red", lty=2) |
|
| 64 |
|
|
| 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", |
|
| 67 |
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) |
|
| 71 |
|
|
| 72 |
plot(lats300, rowMeans(as.matrix(sfd.bg), na.rm=TRUE), type="l", |
|
| 73 |
xlab="Latitude", ylab="Mean flow direction", ylim=ylim) |
|
| 74 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend") |
|
| 75 |
abline(v=60, col="red", lty=2) |
|
| 76 |
|
|
| 77 |
|
|
| 78 |
# |
|
| 79 |
# plot latitudinal profiles of RMSE |
|
| 80 |
# |
|
| 81 |
|
|
| 82 |
# simple helper function to calculate row-wise RMSEs, accounting for the |
|
| 83 |
# fact that flow dir values are circular (0-360), so the difference |
|
| 84 |
# between e.g. 5 and 355 should only be 10 |
|
| 85 |
rmse <- function(r1, r2, na.rm=TRUE, use) {
|
|
| 86 |
diffs <- abs(as.matrix(r1) - as.matrix(r2)) |
|
| 87 |
if (!missing(use)) diffs[!use] <- NA |
|
| 88 |
diffs[] <- ifelse(diffs>180, 360-diffs, diffs) |
|
| 89 |
sqrt(rowMeans(diffs^2, na.rm=na.rm)) |
|
| 90 |
} |
|
| 91 |
|
|
| 92 |
par(mfrow=c(2,3), omi=c(1,1,1,1)) |
|
| 93 |
|
|
| 94 |
ylim <- c(0, 100) |
|
| 95 |
|
|
| 96 |
# ...with respect to ASTER |
|
| 97 |
plot(lats300, rmse(sfd.uncor, sfd.aster), type="l", xlab="Latitude", |
|
| 98 |
ylab="RMSE", ylim=ylim) |
|
| 99 |
lines(lats150, rmse(crop(sfd.uncor, extent(sfd.srtm)), sfd.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), max(ylim)-5, pos=4, font=3, labels="uncorrected") |
|
| 103 |
abline(v=60, col="red", lty=2) |
|
| 104 |
mtext(expression(paste( |
|
| 105 |
"Flowdir discrepancies with respect to separate ASTER/SRTM components (",
|
|
| 106 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 107 |
|
|
| 108 |
#plot(lats300, rmse(sfd.eramp, sfd.aster), type="l", xlab="Latitude", |
|
| 109 |
plot(lats300, rep(0, 300), type="n", xlab="Latitude", |
|
| 110 |
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") |
|
| 114 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="exponential ramp") |
|
| 115 |
text(mean(lats300), mean(ylim), font=3, labels="(skipped)") |
|
| 116 |
#abline(v=60, col="red", lty=2) |
|
| 117 |
|
|
| 118 |
plot(lats300, rmse(sfd.bg, sfd.aster), type="l", xlab="Latitude", |
|
| 119 |
ylab="RMSE", ylim=ylim) |
|
| 120 |
lines(lats150, rmse(crop(sfd.bg, extent(sfd.srtm)), sfd.srtm), col="blue") |
|
| 121 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
|
| 122 |
lty=c(1, 1), bty="n") |
|
| 123 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend") |
|
| 124 |
abline(v=60, col="red", lty=2) |
|
| 125 |
|
|
| 126 |
# ...with respect to CDEM |
|
| 127 |
plot(lats300, rmse(sfd.uncor, sfd.can), type="l", xlab="Latitude", |
|
| 128 |
ylab="RMSE", ylim=ylim) |
|
| 129 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="uncorrected") |
|
| 130 |
abline(v=60, col="red", lty=2) |
|
| 131 |
mtext(expression(paste( |
|
| 132 |
"Flowdir discrepancies with respect to Canada DEM (",
|
|
| 133 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
|
| 134 |
|
|
| 135 |
#plot(lats300, rmse(sfd.eramp, sfd.can), type="l", xlab="Latitude", |
|
| 136 |
plot(lats300, rep(0, 300), type="n", xlab="Latitude", |
|
| 137 |
ylab="RMSE", ylim=ylim) |
|
| 138 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="exponential ramp") |
|
| 139 |
text(mean(lats300), mean(ylim), font=3, labels="(skipped)") |
|
| 140 |
#abline(v=60, col="red", lty=2) |
|
| 141 |
|
|
| 142 |
plot(lats300, rmse(sfd.bg, sfd.can), type="l", xlab="Latitude", |
|
| 143 |
ylab="RMSE", ylim=ylim) |
|
| 144 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend") |
|
| 145 |
abline(v=60, col="red", lty=2) |
|
| 146 |
|
|
| 147 |
# close pdf device driver |
|
| 148 |
dev.off() |
|
| terrain/slope/slope-assessment.R | ||
|---|---|---|
| 61 | 61 |
# |
| 62 | 62 |
|
| 63 | 63 |
# simple helper function to calculate row-wise RMSEs |
| 64 |
rmse <- function(r1, r2, na.rm=TRUE) {
|
|
| 65 |
sqrt(rowMeans(as.matrix((r1 - r2)^2), na.rm=na.rm)) |
|
| 64 |
rmse <- function(r1, r2, na.rm=TRUE, use) {
|
|
| 65 |
diffs <- abs(as.matrix(r1) - as.matrix(r2)) |
|
| 66 |
if (!missing(use)) diffs[!use] <- NA |
|
| 67 |
sqrt(rowMeans(diffs^2, na.rm=na.rm)) |
|
| 66 | 68 |
} |
| 67 | 69 |
|
| 68 | 70 |
par(mfrow=c(2,3), omi=c(1,1,1,1)) |
Also available in: Unified diff
added R code to create various boundary assessment figures