1
|
# R script for ploting a pattern of ASTER2-SRTM deltas as a function of ASTER2 elevation
|
2
|
# The results (plots) are saved at
|
3
|
# "/data/project/organisms/DEM?Yuni/documents/meanElv"
|
4
|
#
|
5
|
# Original author: Jim Regetz
|
6
|
# [13-JUL-2011]
|
7
|
# Edits by Yuina to auto-apply for nine regions.
|
8
|
#
|
9
|
# 9-Dec-2011
|
10
|
# Yuina Nunokawa
|
11
|
|
12
|
|
13
|
|
14
|
|
15
|
|
16
|
|
17
|
# plot pattern of ASTER2-SRTM deltas as a function of ASTER2 elevation
|
18
|
plotMeanDeltaByElev <- function(delta.vals, elev, ...) {
|
19
|
mean.by.elev <- tapply(delta.vals, elev, mean)
|
20
|
sd.by.elev <- tapply(delta.vals, elev, sd)
|
21
|
n.by.elev <- tapply(delta.vals, elev, length)
|
22
|
se.by.elev <- sd.by.elev/sqrt(n.by.elev)
|
23
|
na.se.points <- mean.by.elev[is.na(se.by.elev)]
|
24
|
se.by.elev[is.na(se.by.elev)] <- 0
|
25
|
elev <- as.numeric(names(mean.by.elev))
|
26
|
plot(elev, mean.by.elev, pch=16, xlim=c(0, max(elev)),
|
27
|
ylim=c(min(mean.by.elev - se.by.elev), max(mean.by.elev + se.by.elev)),
|
28
|
type="n", ...)
|
29
|
segments(elev, mean.by.elev-se.by.elev,
|
30
|
as.numeric(names(mean.by.elev)), mean.by.elev+se.by.elev,col="grey")
|
31
|
points(elev, mean.by.elev, pch=".")
|
32
|
points(as.numeric(names(na.se.points)), na.se.points, pch=4,
|
33
|
col="red", cex=0.5)
|
34
|
}
|
35
|
|
36
|
|
37
|
|
38
|
|
39
|
plotDeltaBins <- function(delta.vals, elev, bin.min, bin.width, bin.max=1500,
|
40
|
outline=FALSE, ...) {
|
41
|
breaks <- seq(bin.min, bin.max, by=bin.width)
|
42
|
midpts <- c(paste("<", bin.min, sep=""), head(breaks, -1) + bin.width/2,
|
43
|
paste(">", bin.max, sep=""))
|
44
|
elev <- cut(elev, breaks=c(0, breaks, Inf), labels=midpts)
|
45
|
bp <- boxplot(delta.vals ~ elev, outline=outline, col="lightgray",
|
46
|
frame=FALSE, ...)
|
47
|
text(1:length(bp$n), bp$stats[5,], labels=round(bp$n/1000),
|
48
|
pos=3, cex=0.5, offset=0.2, font=3, col="gray")
|
49
|
#axis(3, at=seq_along(bp$n), labels=paste(round(bp$n/1000), "k", sep=""),
|
50
|
# cex.axis=0.7, tick=FALSE, font=3, line=-1)
|
51
|
#mtext("n =", side=3, adj=0, font=3, cex=0.7)
|
52
|
abline(h=median(delta.vals), col="red", lty=2)
|
53
|
invisible(bp)
|
54
|
}
|
55
|
|
56
|
|
57
|
|
58
|
|
59
|
library(raster)
|
60
|
|
61
|
aster.dir <- "/data/project/organisms/DEM/Yuni/Data/aster2"
|
62
|
srtm.dir <- "/data/project/organisms/DEM/Yuni/Data/srtm"
|
63
|
aster.straddle <- list.files(aster.dir, pattern="^aster2.*_straddle.tif")
|
64
|
srtm.below <- list.files(srtm.dir, pattern="^srtm_.*_below.tif")
|
65
|
fused.straddle <- list.files(aster.dir, pattern="^fused.*_straddle.tif")
|
66
|
fused.bgau <- list.files(aster.dir, pattern="^fused.*_blendgau.tif")
|
67
|
|
68
|
|
69
|
|
70
|
|
71
|
|
72
|
deltas <- function(d.aster, d.srtm, d.uncor, d.bg, file.name, aster.dir, srtm.dir){
|
73
|
|
74
|
# initialize output pdf device driver
|
75
|
name <- paste("deltas_", file.name, ".pdf", sep="")
|
76
|
pdf(file.path("/data/project/organisms/DEM/Yuni/documents/deltas", file=name), height=8, width=11.5)
|
77
|
|
78
|
|
79
|
# plot mean delta elevations
|
80
|
d.aster2.crop.vals <- values(crop(d.aster, extent(d.srtm)))
|
81
|
d.srtm.vals <- values(d.srtm)
|
82
|
delta.vals <- d.aster2.crop.vals - d.srtm.vals
|
83
|
plotMeanDeltaByElev(delta.vals, d.aster2.crop.vals,
|
84
|
xlab="ASTER2 elevation (m)", ylab="ASTER2-SRTM difference (m)")
|
85
|
|
86
|
# plot aster2 srtm bins
|
87
|
d.aster2.crop.vals <- values(crop(d.aster, extent(d.srtm)))
|
88
|
d.srtm.vals <- values(d.srtm)
|
89
|
delta.vals <- d.aster2.crop.vals - d.srtm.vals
|
90
|
# d.aster2.crop.vals <- d.aster2.crop.vals[d.srtm.vals>0]
|
91
|
# d.srtm.vals <- d.srtm.vals[d.srtm.vals>0]
|
92
|
|
93
|
plotDeltaBins(delta.vals, d.srtm.vals, 150, 50, 1500, las=2,
|
94
|
cex.axis=0.8, xlab="Midpoints of SRTM elevation bins (m)",
|
95
|
ylab="ASTER2 - SRTM difference (m)")
|
96
|
|
97
|
|
98
|
# plot scatter of aster vs srtm
|
99
|
plot(jitter(d.srtm.vals), jitter(d.aster2.crop.vals), pch=".",
|
100
|
xlab="SRTM elevation (m)", ylab="ASTER2 elevation (m)", cex=0.5)
|
101
|
abline(median(delta.vals), 1, col="red", cex=0.5)
|
102
|
abline(0, 1, col="blue", lty=2, cex=0.5)
|
103
|
|
104
|
# add inset histogram of differences
|
105
|
opar <- par(fig=c(0.55, 0.95, 0.1, 0.6), new=TRUE)
|
106
|
h <- hist(delta.vals[abs(delta.vals)<60], breaks=48, xlab=NA, main=NULL,
|
107
|
col=grey(0.8), border=grey(0.3), yaxt="n", ylab=NA, cex.axis=0.5,
|
108
|
cex.lab=0.5, tcl=-0.25, mgp=c(3,0,0))
|
109
|
text(10, 0.4*max(h$counts), labels=paste("Entire range:\n(",
|
110
|
min(delta.vals), ", ", max(delta.vals), ")", sep=""), cex=0.6, adj=c(0,0))
|
111
|
mtext("ASTER2 - SRTM (m)", side=1, cex=0.5, line=0.6)
|
112
|
|
113
|
dev.off ()
|
114
|
}
|
115
|
|
116
|
|
117
|
|
118
|
num.regions <- 9
|
119
|
plot.deltas <-lapply(1:num.regions, function(region, aster.dir, srtm.dir){
|
120
|
d.aster <- raster(file.path(aster.dir, aster.straddle[region]))
|
121
|
d.srtm <- raster(file.path(srtm.dir, srtm.below[region]))
|
122
|
d.uncor <- raster(file.path(aster.dir, fused.straddle[region]))
|
123
|
d.bg <- raster(file.path(aster.dir, fused.bgau[region]))
|
124
|
file.name <- substr(srtm.below[region],6, 13)
|
125
|
|
126
|
deltas(d.aster, d.srtm, d.uncor, d.bg, file.name, aster.dir, srtm.dir)
|
127
|
}, aster.dir=aster.dir, srtm.dir=srtm.dir)
|
128
|
|
129
|
|