Revision 7cf747b0
Added by Jim Regetz over 13 years ago
- ID 7cf747b082743603bbcf9022e5ed1303b448bb4b
| terrain/aspect/aspect-assessment.R | ||
|---|---|---|
| 11 | 11 |
datadir <- "/home/regetz/media/temp/terrain/aspect" |
| 12 | 12 |
|
| 13 | 13 |
# load aspect rasters |
| 14 |
s.aster <- raster(file.path(datadir, "aster_300straddle_a.tif"))
|
|
| 15 |
s.srtm <- raster(file.path(datadir, "srtm_150below_a.tif"))
|
|
| 16 |
s.uncor <- raster(file.path(datadir, "fused_300straddle_a.tif"))
|
|
| 17 |
s.eramp <- raster(file.path(datadir, "fused_300straddle_rampexp_a.tif"))
|
|
| 18 |
s.bg <- raster(file.path(datadir, "fused_300straddle_blendgau_a.tif"))
|
|
| 19 |
s.can <- raster(file.path(datadir, "cdem_300straddle_a.tif"))
|
|
| 14 |
a.aster <- raster(file.path(datadir, "aster_300straddle_a.tif"))
|
|
| 15 |
a.srtm <- raster(file.path(datadir, "srtm_150below_a.tif"))
|
|
| 16 |
a.uncor <- raster(file.path(datadir, "fused_300straddle_a.tif"))
|
|
| 17 |
a.enblend <- raster(file.path(datadir, "fused_300straddle_enblend_a.tif"))
|
|
| 18 |
a.bg <- raster(file.path(datadir, "fused_300straddle_blendgau_a.tif"))
|
|
| 19 |
a.can <- raster(file.path(datadir, "cdem_300straddle_a.tif"))
|
|
| 20 | 20 |
|
| 21 | 21 |
# extract raster latitudes for later |
| 22 |
lats300 <- yFromRow(s.aster, 1:nrow(s.aster))
|
|
| 23 |
lats150 <- yFromRow(s.srtm, 1:nrow(s.srtm))
|
|
| 22 |
lats300 <- yFromRow(a.aster, 1:nrow(a.aster))
|
|
| 23 |
lats150 <- yFromRow(a.srtm, 1:nrow(a.srtm))
|
|
| 24 | 24 |
|
| 25 | 25 |
# initialize output pdf device driver |
| 26 | 26 |
pdf("aspect-assessment.pdf", height=8, width=11.5)
|
| ... | ... | |
| 34 | 34 |
|
| 35 | 35 |
ylim <- c(160, 180) |
| 36 | 36 |
|
| 37 |
plot(lats300, rowMeans(as.matrix(s.uncor), na.rm=TRUE), type="l",
|
|
| 37 |
plot(lats300, rowMeans(as.matrix(a.can), na.rm=TRUE), type="l",
|
|
| 38 | 38 |
xlab="Latitude", ylab="Mean aspect", ylim=ylim) |
| 39 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="uncorrected") |
|
| 39 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Original DEMs") |
|
| 40 |
lines(lats300, rowMeans(as.matrix(a.aster), na.rm=TRUE), col="blue") |
|
| 41 |
lines(lats150, rowMeans(as.matrix(a.srtm), na.rm=TRUE), col="red") |
|
| 42 |
legend("bottomright", legend=c("ASTER", "SRTM", "CDED"), col=c("blue",
|
|
| 43 |
"red", "black"), lty=c(1, 1), bty="n") |
|
| 40 | 44 |
abline(v=60, col="red", lty=2) |
| 41 | 45 |
mtext(expression(paste("Latitudinal profiles of mean aspect (",
|
| 42 | 46 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
| 43 | 47 |
|
| 44 |
plot(lats300, rowMeans(as.matrix(s.can), na.rm=TRUE), type="l",
|
|
| 48 |
plot(lats300, rowMeans(as.matrix(a.uncor), na.rm=TRUE), type="l",
|
|
| 45 | 49 |
xlab="Latitude", ylab="Mean aspect", ylim=ylim) |
| 46 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="Canada DEM")
|
|
| 50 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="simple fused")
|
|
| 47 | 51 |
abline(v=60, col="red", lty=2) |
| 48 | 52 |
|
| 49 |
plot(lats300, rowMeans(as.matrix(s.eramp), na.rm=TRUE), type="l",
|
|
| 53 |
plot(lats300, rowMeans(as.matrix(a.enblend), na.rm=TRUE), type="l",
|
|
| 50 | 54 |
xlab="Latitude", ylab="Mean aspect", ylim=ylim) |
| 51 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="exponential ramp")
|
|
| 55 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="multires spline")
|
|
| 52 | 56 |
abline(v=60, col="red", lty=2) |
| 53 | 57 |
|
| 54 |
plot(lats300, rowMeans(as.matrix(s.bg), na.rm=TRUE), type="l",
|
|
| 58 |
plot(lats300, rowMeans(as.matrix(a.bg), na.rm=TRUE), type="l",
|
|
| 55 | 59 |
xlab="Latitude", ylab="Mean aspect", ylim=ylim) |
| 56 | 60 |
text(min(lats300), min(ylim)+0.5, pos=4, font=3, labels="gaussian blend") |
| 57 | 61 |
abline(v=60, col="red", lty=2) |
| ... | ... | |
| 76 | 80 |
ylim <- c(0, 100) |
| 77 | 81 |
|
| 78 | 82 |
# ...with respect to ASTER |
| 79 |
plot(lats300, rmse(s.uncor, s.aster), type="l", xlab="Latitude",
|
|
| 83 |
plot(lats300, rmse(a.uncor, a.aster), type="l", xlab="Latitude",
|
|
| 80 | 84 |
ylab="RMSE", ylim=ylim) |
| 81 |
lines(lats150, rmse(crop(s.uncor, extent(s.srtm)), s.srtm), col="blue")
|
|
| 85 |
lines(lats150, rmse(crop(a.uncor, extent(a.srtm)), a.srtm), col="blue")
|
|
| 82 | 86 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
| 83 | 87 |
lty=c(1, 1), bty="n") |
| 84 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="uncorrected")
|
|
| 88 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="simple fused")
|
|
| 85 | 89 |
abline(v=60, col="red", lty=2) |
| 86 | 90 |
mtext(expression(paste( |
| 87 | 91 |
"Aspect discrepancies with respect to separate ASTER/SRTM components (",
|
| 88 | 92 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
| 89 | 93 |
|
| 90 |
plot(lats300, rmse(s.eramp, s.aster), type="l", xlab="Latitude",
|
|
| 94 |
plot(lats300, rmse(a.enblend, a.aster), type="l", xlab="Latitude",
|
|
| 91 | 95 |
ylab="RMSE", ylim=ylim) |
| 92 |
lines(lats150, rmse(crop(s.eramp, extent(s.srtm)), s.srtm), col="blue")
|
|
| 96 |
lines(lats150, rmse(crop(a.enblend, extent(a.srtm)), a.srtm), col="blue")
|
|
| 93 | 97 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
| 94 | 98 |
lty=c(1, 1), bty="n") |
| 95 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="exponential ramp")
|
|
| 99 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="multires spline")
|
|
| 96 | 100 |
abline(v=60, col="red", lty=2) |
| 97 | 101 |
|
| 98 |
plot(lats300, rmse(s.bg, s.aster), type="l", xlab="Latitude",
|
|
| 102 |
plot(lats300, rmse(a.bg, a.aster), type="l", xlab="Latitude",
|
|
| 99 | 103 |
ylab="RMSE", ylim=ylim) |
| 100 |
lines(lats150, rmse(crop(s.bg, extent(s.srtm)), s.srtm), col="blue")
|
|
| 104 |
lines(lats150, rmse(crop(a.bg, extent(a.srtm)), a.srtm), col="blue")
|
|
| 101 | 105 |
legend("topright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
| 102 | 106 |
lty=c(1, 1), bty="n") |
| 103 | 107 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend") |
| 104 | 108 |
abline(v=60, col="red", lty=2) |
| 105 | 109 |
|
| 106 | 110 |
# ...with respect to CDEM |
| 107 |
plot(lats300, rmse(s.uncor, s.can), type="l", xlab="Latitude",
|
|
| 111 |
plot(lats300, rmse(a.uncor, a.can), type="l", xlab="Latitude",
|
|
| 108 | 112 |
ylab="RMSE", ylim=ylim) |
| 109 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="uncorrected")
|
|
| 113 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="simple fused")
|
|
| 110 | 114 |
abline(v=60, col="red", lty=2) |
| 111 | 115 |
mtext(expression(paste( |
| 112 | 116 |
"Aspect discrepancies with respect to Canada DEM (",
|
| 113 | 117 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
| 114 | 118 |
|
| 115 |
plot(lats300, rmse(s.eramp, s.can), type="l", xlab="Latitude",
|
|
| 119 |
plot(lats300, rmse(a.enblend, a.can), type="l", xlab="Latitude",
|
|
| 116 | 120 |
ylab="RMSE", ylim=ylim) |
| 117 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="exponential ramp")
|
|
| 121 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="multires spline")
|
|
| 118 | 122 |
abline(v=60, col="red", lty=2) |
| 119 | 123 |
|
| 120 |
plot(lats300, rmse(s.bg, s.can), type="l", xlab="Latitude",
|
|
| 124 |
plot(lats300, rmse(a.bg, a.can), type="l", xlab="Latitude",
|
|
| 121 | 125 |
ylab="RMSE", ylim=ylim) |
| 122 | 126 |
text(min(lats300), max(ylim)-5, pos=4, font=3, labels="gaussian blend") |
| 123 | 127 |
abline(v=60, col="red", lty=2) |
| ... | ... | |
| 147 | 151 |
ylim <- c(0, 1) |
| 148 | 152 |
|
| 149 | 153 |
# ...with respect to ASTER |
| 150 |
plot(lats300, corByLat(s.uncor, s.aster), type="l", xlab="Latitude",
|
|
| 154 |
plot(lats300, corByLat(a.uncor, a.aster), type="l", xlab="Latitude",
|
|
| 151 | 155 |
ylab="Correlation", ylim=ylim) |
| 152 |
lines(lats150, corByLat(crop(s.uncor, extent(s.srtm)), s.srtm), col="blue")
|
|
| 156 |
lines(lats150, corByLat(crop(a.uncor, extent(a.srtm)), a.srtm), col="blue")
|
|
| 153 | 157 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
| 154 | 158 |
lty=c(1, 1), bty="n") |
| 155 |
text(min(lats300), min(ylim), pos=4, font=3, labels="uncorrected")
|
|
| 159 |
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused")
|
|
| 156 | 160 |
abline(v=60, col="red", lty=2) |
| 157 | 161 |
mtext(expression(paste( |
| 158 | 162 |
"Aspect correlations with respect to separate ASTER/SRTM components (",
|
| 159 | 163 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
| 160 | 164 |
|
| 161 |
plot(lats300, corByLat(s.eramp, s.aster), type="l", xlab="Latitude",
|
|
| 165 |
plot(lats300, corByLat(a.enblend, a.aster), type="l", xlab="Latitude",
|
|
| 162 | 166 |
ylab="Correlation", ylim=ylim) |
| 163 |
lines(lats150, corByLat(crop(s.eramp, extent(s.srtm)), s.srtm), col="blue")
|
|
| 167 |
lines(lats150, corByLat(crop(a.enblend, extent(a.srtm)), a.srtm), col="blue")
|
|
| 164 | 168 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
| 165 | 169 |
lty=c(1, 1), bty="n") |
| 166 |
text(min(lats300), min(ylim), pos=4, font=3, labels="exponential ramp")
|
|
| 170 |
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")
|
|
| 167 | 171 |
abline(v=60, col="red", lty=2) |
| 168 | 172 |
|
| 169 |
plot(lats300, corByLat(s.bg, s.aster), type="l", xlab="Latitude",
|
|
| 173 |
plot(lats300, corByLat(a.bg, a.aster), type="l", xlab="Latitude",
|
|
| 170 | 174 |
ylab="Correlation", ylim=ylim) |
| 171 |
lines(lats150, corByLat(crop(s.bg, extent(s.srtm)), s.srtm), col="blue")
|
|
| 175 |
lines(lats150, corByLat(crop(a.bg, extent(a.srtm)), a.srtm), col="blue")
|
|
| 172 | 176 |
legend("bottomright", legend=c("ASTER", "SRTM"), col=c("black", "blue"),
|
| 173 | 177 |
lty=c(1, 1), bty="n") |
| 174 | 178 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
| 175 | 179 |
abline(v=60, col="red", lty=2) |
| 176 | 180 |
|
| 177 | 181 |
# ...with respect to CDEM |
| 178 |
plot(lats300, corByLat(s.uncor, s.can), type="l", xlab="Latitude",
|
|
| 182 |
plot(lats300, corByLat(a.uncor, a.can), type="l", xlab="Latitude",
|
|
| 179 | 183 |
ylab="Correlation", ylim=ylim) |
| 180 |
text(min(lats300), min(ylim), pos=4, font=3, labels="uncorrected")
|
|
| 184 |
text(min(lats300), min(ylim), pos=4, font=3, labels="simple fused")
|
|
| 181 | 185 |
abline(v=60, col="red", lty=2) |
| 182 | 186 |
mtext(expression(paste( |
| 183 | 187 |
"Aspect correlations with respect to Canada DEM (",
|
| 184 | 188 |
136*degree, "W to ", 96*degree, "W)")), adj=0, line=2, font=2) |
| 185 | 189 |
|
| 186 |
plot(lats300, corByLat(s.eramp, s.can), type="l", xlab="Latitude",
|
|
| 190 |
plot(lats300, corByLat(a.enblend, a.can), type="l", xlab="Latitude",
|
|
| 187 | 191 |
ylab="Correlation", ylim=ylim) |
| 188 |
text(min(lats300), min(ylim), pos=4, font=3, labels="exponential ramp")
|
|
| 192 |
text(min(lats300), min(ylim), pos=4, font=3, labels="multires spline")
|
|
| 189 | 193 |
abline(v=60, col="red", lty=2) |
| 190 | 194 |
|
| 191 |
plot(lats300, corByLat(s.bg, s.can), type="l", xlab="Latitude",
|
|
| 195 |
plot(lats300, corByLat(a.bg, a.can), type="l", xlab="Latitude",
|
|
| 192 | 196 |
ylab="Correlation", ylim=ylim) |
| 193 | 197 |
text(min(lats300), min(ylim), pos=4, font=3, labels="gaussian blend") |
| 194 | 198 |
abline(v=60, col="red", lty=2) |
Also available in: Unified diff
updated boundary assessment code to use multires spline results