1 |
7526fb1c
|
Jim Regetz
|
# 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 |
|
|
|