Project

General

Profile

Download (8.62 KB) Statistics
| Branch: | Revision:
1
##############################################################################################
2
#
3
# CreateImageDiffPlots
4
# 
5
# R script generates collection plots that explore elevation discontinuities 
6
# border between ASTER and STRM data in mosaic images under development for Iplant.
7
# Plots overwritten to single graphics device for first version. 
8
# Panel Plots to be added to next version.
9
#
10
# To run:
11
#
12
#  1) place this file and the input file "tableForMark4000_5_8_SortColID.csv" in folder
13
#  2) start R
14
#  3) > source("MosaicCDEMDifferencePlots.r")
15
#
16
# Input file: Csv table containing elevations for randomly-selected point pairs from 
17
# ASTER/CGIAR mosaic and CDEM images. Table is generated by the script: 'makeImagePairTable.r'
18
# Author: Rick Reeves, NCEAS
19
# May 1, 2011
20
##############################################################################################
21
#CreateImageDiffPlots <- function()
22
#{
23

    
24
# Read input table that was sorted OFFLINE in Excel on the first column (ColumnID)
25

    
26
   pointTable <-read.csv("tableForMark4000_5_8_SortColID.csv")
27
#
28
# Table created by randomly sampling from two superimposed 
29
# image: 
30
#  1) DOM mosaic image comprised of ASTER and SRTM components.
31
#  2) 'Baseline' Canadian DEM (CDEM) image.
32
# 
33
# The input table contains three rows for each randomly-selected
34
# pixel pair from both images. Each row contains two pixel pairs,
35
# the first pair drawn from the image mosaic, the second pair
36
# drawn from the CDEM image: 
37
#   First pair: North pixel, South pixel (ASTER/SRTM mosaic)
38
#   Second pair: North pixel, South pixel (CDEM)
39
#
40
#  The first row of each 'triplet' contains pixel pairs North of border,
41
#  The second row contains pixel pairs spanning  border,
42
#  The third row contains pixel pairs South of border,
43
#
44
# This script generates a series of plots that display
45
# differences between:
46
#    1) The mosaic and CDEM images
47
#    2) Image pixels on and away from the ASTER / SRTM boundary.
48
# 
49
   northRowIndexes = seq(from=1, to=(nrow(pointTable) - 3),by=3)
50
   borderRowIndexes = seq(from=2, to=(nrow(pointTable) - 1),by=3)   
51
   southRowIndexes = seq(from=3, to=(nrow(pointTable)),by=3)
52
#
53
# calculate and append the difference between elevations
54
# and CDEM
55
# these lines create the inputs for differnce image plots for each of three
56
# pixel pair subsets: North of border (All Aster), border (combo Aster/Srtm),
57
#                     South of border (all Srtm)
58
# 
59
# First, add the 'difference columns' to the entire table
60
#
61
   pointTable <-cbind(pointTable,(pointTable$elevNorth - pointTable$elevSouth))
62
   pointTable <-cbind(pointTable,(pointTable$cdemNorth - pointTable$cdemSouth))
63
   pointTable <-cbind(pointTable,(pointTable$elevNorth - pointTable$cdemNorth))
64
   pointTable <-cbind(pointTable,(pointTable$elevSouth - pointTable$cdemSouth))
65
#   
66
   colnames(pointTable)[6] <- "diffMosaicNorthSouth"
67
   colnames(pointTable)[7] <- "diffCDEMNorthSouth"
68
   colnames(pointTable)[8] <- "diffNorthMosaicCDEM"
69
   colnames(pointTable)[9] <- "diffSouthMosaicCDEM"   
70

    
71
# Difference between Mosaic (ASTER or CGIAR or border) 
72
# and CDEM elevation as pertentage of the mosaic elevation
73

    
74
   pointTable <-cbind(pointTable,(pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100)) 
75
   pointTable <-cbind(pointTable,(pointTable$diffSouthMosaicCDEM/pointTable$elevSouth * 100)) 
76
   
77
   colnames(pointTable)[10] <- "magDiffMosaicCDEMNorthPct"
78
   colnames(pointTable)[11] <- "magDiffMosaicCDEMSouthPct"   
79

    
80
# For the plots, subdivide the table into three segments: 
81
#     rows north of border
82
#     rows crossing border
83
#     rows south of border
84

    
85
   northRowTbl = pointTable[northRowIndexes,]
86
   borderRowTbl = pointTable[borderRowIndexes,]
87
   southRowTbl = pointTable[southRowIndexes,]
88

    
89
#
90
# Create plots
91
#
92
# First, plot North/South pixel elevation for North/Border/South subsets 'distance East from Western edge'.
93
#
94
   plot(northRowTbl$ColumnID,northRowTbl$elevNorth,main="North (ASTER) Mosaic North Pixel Elevation",col="red",pch=".")
95
message("Hit key for next plot...")
96
browser()   
97
   plot(borderRowTbl$ColumnID,borderRowTbl$elevNorth,main="Border Mosaic North Elevation",col="darkgreen",pch=".")
98
message("Hit key for next plot...")
99
browser()   
100
   plot(southRowTbl$ColumnID,southRowTbl$elevNorth,main="South (SRTM) Mosaic North Elevation",col="blue",pch=".")
101
message("Hit key for next plot...")
102
browser()   
103
   plot(northRowTbl$ColumnID,northRowTbl$elevNorth,main="North (ASTER) CDEM North Pixel Elevation",col="red",pch=".")
104
message("Hit key for next plot...")
105
browser()   
106
   plot(borderRowTbl$ColumnID,borderRowTbl$elevNorth,main="Border CDEM North Pixel Elevation",col="darkgreen",pch=".")
107
message("Hit key for next plot...")
108
browser()   
109
   plot(southRowTbl$ColumnID,southRowTbl$elevNorth,main="South (SRTM) CDEM North Elevation",col="blue",pch=".")
110
message("Hit key for next plot...")
111
browser()
112
   plot(northRowTbl$ColumnID,northRowTbl$elevSouth,main="North (ASTER) Mosaic South Pixel Elevation",col="red",pch=".")
113
message("Hit key for next plot...")
114
browser()   
115
   plot(borderRowTbl$ColumnID,borderRowTbl$elevSouth,main="Border Mosaic South Elevation",col="darkgreen",pch=".")
116
message("Hit key for next plot...")
117
browser()   
118
   plot(southRowTbl$ColumnID,southRowTbl$elevSouth,main="South (SRTM) Mosaic South Elevation",col="blue",pch=".")
119
message("Hit key for next plot...")
120
browser()
121
   plot(northRowTbl$ColumnID,northRowTbl$elevSouth,main="North (ASTER) CDEM South Pixel Elevation",col="red",pch=".")
122
message("Hit key for next plot...")
123
browser()   
124
   plot(borderRowTbl$ColumnID,borderRowTbl$elevSouth,main="Border CDEM South Pixel Elevation",col="darkgreen",pch=".")
125
message("Hit key for next plot...")
126
browser()   
127
   plot(southRowTbl$ColumnID,southRowTbl$elevSouth,main="South (SRTM) CDEM South Elevation",col="blue",pch=".")   
128
message("Hit key for next plot...")
129
browser()   
130
#
131
# Next, (for now) difference between North and South pixels in pair, for various regions
132
# 
133
   plot(northRowTbl$ColumnID,northRowTbl$diffMosaicNorthSouth,main="North (ASTER) Mosaic North/South Diff",col="red",pch=".")
134
message("Hit key for next plot...")
135
browser()   
136
   plot(borderRowTbl$ColumnID,borderRowTbl$diffMosaicNorthSouth,main="Border Mosaic North/South Diff",col="darkgreen",pch=".")
137
message("Hit key for next plot...")
138
browser()   
139
   plot(southRowTbl$ColumnID,southRowTbl$diffMosaicNorthSouth,main="South (SRTM) Mosaic North/South Diff",col="darkgreen",pch=".")   
140
message("Hit key for next plot...")
141
browser()   
142
#   
143
   plot(northRowTbl$ColumnID,northRowTbl$diffCDEMNorthSouth,main="North (ASTER) CDEM North/South Diff",col="red",pch=".")
144
message("Hit key for next plot...")
145
browser()   
146
   plot(borderRowTbl$ColumnID,borderRowTbl$diffCDEMNorthSouth,main="Border CDEM North/South Diff",col="darkgreen",pch=".")
147
message("Hit key for next plot...")
148
browser()   
149
   plot(southRowTbl$ColumnID,southRowTbl$diffCDEMNorthSouth,main="South (SRTM) CDEM North/South Diff",col="blue",pch=".")
150
message("Hit key for next plot...")
151
browser()   
152

    
153
# Next: the Mosaic/CDEM difference divided by the (ASTER or SRTM) elevation
154
  
155
  plot(northRowTbl$ColumnID,northRowTbl$magDiffMosaicCDEMNorthPct,col="red",pch=".",main="Elev Diff as PerCent of Mosaic Elev (North)")
156
message("Hit key for next plot...")
157
browser()  
158
  plot(borderRowTbl$ColumnID,borderRowTbl$magDiffMosaicCDEMNorthPct,col="darkgreen",pch=".",main="Elev Diff as PerCent of Mosaic Elev (Border)")  
159
message("Hit key for next plot...")
160
browser()  
161
  plot(southRowTbl$ColumnID,southRowTbl$magDiffMosaicCDEMSouthPct,col="blue",pch=".",main="Elev Diff as PerCent of Mosaic Elev (South)")  
162
message("Hit key for next plot...")
163
browser()  
164
# Finally: Mosaic vs CDEM difference vs elevation (north or south)
165

    
166
  plot(northRowTbl$elevNorth,northRowTbl$diffCDEMNorthSouth,col="red",pch=".",main="ASTER Elev Diff vs Mosaic Elev (North Pixel)")
167
# plot(northRowTbl$elevSouth,southRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="ASTER Elev Diff vs Mosaic Elev (South Pixel)")  
168
message("Hit key for next plot...")
169
browser()  
170
  plot(borderRowTbl$elevNorth,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="Border Elev Diff vs Mosaic Elev (North Pixel)")
171
#  plot(borderRowTbl$elevSouth,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="Border Elev Diff vs Mosaic Elev (South Pixel)")  
172
message("Hit key for next plot...")
173
browser()  
174
  plot(southRowTbl$elevNorth,southRowTbl$diffCDEMNorthSouth,col="blue",pch=".",main="SRTM Elev Diff vs Mosaic Elev (North Pixel)")
175
#  plot(southRowTbl$elevSouth,southRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="SRTM Elev Diff vs Mosaic Elev (South Pixel)")  
176
message("end of plots")
177
#}
(8-8/19)