Project

General

Profile

Download (16.2 KB) Statistics
| Branch: | Revision:
1
##############################################################################################
2
#
3
# CreateImageDiffPlots
4
# 
5
# R script generates collection of nine plots that display the distributions of  populations of 
6
# three elevation pixel pairs in proximity to the border between ASTER and STRM data in mosaic 
7
# images# under development for the Environment and Organisms project.
8
# 
9
#
10
# Inputs: 
11
#   1) Comma Separated Value (CSV) table containing randomly-sampled elevation value pairs,
12
#      extracted from ASTER/CGIAR mosaic and CDEM images, created by the R script: makeImagePairTable.r
13
#      Note: At present, be sure that this file has been sorted by column 1 (ColumnID) in Excel
14
#      before using in this program. 
15
#      NOTE: This filename is set in the body of the script.
16
#
17
#   2) sampling factor: integer (range=1 to number of recors in input table. Determines
18
#      the size of randomly-selected sampe of total table record 'triplets' (north, border,
19
#      and south) rows of each column sample) to be displayed in the plots: 
20
#           sampling factor = 1  : plot every table record
21
#                             10 : plot a random sample containing 1/10th of records
22
#                            100 : plot random sample containing 1/100th of records.
23
#   3) JpgPlotFileFlag: Set to TRUE to have *individual* plots written to JPG files.
24
#
25
# To run:
26
#
27
#  1) place this file and the input file "tableForMark4000_5_8_SortColID.csv" in folder
28
#  2) start R
29
#  3) > source("CreateImageDiffPlotsOverlay.r")
30
#  4) > CreateImageDiffPlots(sampling factor)
31
#     < press <cr> key to view plots sequentially
32
#  
33
# TODO: add symbol legend to each plot.
34
# Author: Rick Reeves, NCEAS
35
# May 14, 2011
36
# May 17, 2011: This version generates 'delta scatterplots' specified by Mark and Jim.
37
# May 18, 2011: If JpgPlotFlag set, write individual plot files to JPG files.
38
##############################################################################################
39
CreateImageDiffPlots <- function(plotSampFact = 1,JpgPlotFileFlag = TRUE)
40
{
41

    
42
# Check plotSampleFact range
43

    
44
   if (plotSampFact < 1)
45
      plotSampFact = 1
46
      
47
# Read input table that was sorted OFFLINE in Excel on the first column (ColumnID)
48

    
49
   pointTable <-read.csv("pixelPairs36000_5_8_TS.csv")
50
# Table created by randomly sampling from two superimposed 
51
# image: 
52
#  1) DOM mosaic image comprised of ASTER and SRTM components.
53
#  2) 'Baseline' Canadian DEM (CDEM) image.
54
# 
55
# The input table contains three rows for each randomly-selected
56
# pixel pair from both images. Each row contains two pixel pairs,
57
# the first pair drawn from the image mosaic, the second pair
58
# drawn from the CDEM image: 
59
#   First pair: North pixel, South pixel (ASTER/SRTM mosaic)
60
#   Second pair: North pixel, South pixel (CDEM)
61
#
62
#  The first row of each 'triplet' contains pixel pairs North of border,
63
#  The second row contains pixel pairs spanning  border,
64
#  The third row contains pixel pairs South of border,
65
#
66
# This script generates a series of plots that display
67
# differences between:
68
#    1) The mosaic and CDEM images
69
#    2) Image pixels on and away from the ASTER / SRTM boundary.
70
# 
71
   northRowIndexes = seq(from=1, to=(nrow(pointTable) - 3),by=3)
72
   borderRowIndexes = seq(from=2, to=(nrow(pointTable) - 1),by=3)   
73
   southRowIndexes = seq(from=3, to=(nrow(pointTable)),by=3)
74
#
75
# calculate and append the difference between elevations
76
# and CDEM
77
# these lines create the inputs for differnce image plots for each of three
78
# pixel pair subsets: North of border (All Aster), border (combo Aster/Srtm),
79
#                     South of border (all Srtm)
80
# 
81
# First, add the 'difference columns' to the entire table
82
#
83
   pointTable <-cbind(pointTable,(pointTable$elevNorth - pointTable$elevSouth))
84
   pointTable <-cbind(pointTable,(pointTable$cdemNorth - pointTable$cdemSouth))
85
   pointTable <-cbind(pointTable,(pointTable$elevNorth - pointTable$cdemNorth))
86
   pointTable <-cbind(pointTable,(pointTable$elevSouth - pointTable$cdemSouth))
87
#   
88
   colnames(pointTable)[6] <- "diffMosaicNorthSouth"
89
   colnames(pointTable)[7] <- "diffCDEMNorthSouth"
90
   colnames(pointTable)[8] <- "diffNorthMosaicCDEM"
91
   colnames(pointTable)[9] <- "diffSouthMosaicCDEM"   
92
   
93
# add a placeholder for the 'boundary' value, across each table
94

    
95
#  pointTable <-cbind(pointTable,(-1))
96
   
97
#  colnames(pointTable)[10] <- "deltaAcrossBorder"
98

    
99
# Difference between Mosaic (ASTER or CGIAR or border) 
100
# and CDEM elevation as pertentage of the mosaic elevation
101

    
102
   pointTable <-cbind(pointTable,(pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100)) 
103
   pointTable <-cbind(pointTable,(pointTable$diffSouthMosaicCDEM/pointTable$elevSouth * 100)) 
104
   
105
   colnames(pointTable)[10] <- "magDiffMosaicCDEMNorthPct"
106
   colnames(pointTable)[11] <- "magDiffMosaicCDEMSouthPct"   
107

    
108
# For the plots, subdivide the table into three segments: 
109
#     rows north of border
110
#     rows crossing border
111
#     rows south of border
112

    
113
   northRowTblAll = pointTable[northRowIndexes,]
114
   borderRowTblAll = pointTable[borderRowIndexes,]
115
   southRowTblAll = pointTable[southRowIndexes,]
116

    
117
   subset <- 1:nrow(northRowTblAll)
118
   
119
   randSub <- sample(subset,as.integer(length(subset) / plotSampFact)) 
120
   
121
   northRowTbl <- northRowTblAll[randSub,]
122
   borderRowTbl <- borderRowTblAll[randSub,]
123
   southRowTbl <- southRowTblAll[randSub,]   
124
   
125
message("hit key to create each plot...")
126
browser()
127

    
128
# Three plotting characters used
129

    
130
   plotCh1 <- 17 # 'north' (aster) points
131
   plotCh2 <- 18 # 'border' (aster+srtm) points 
132
   plotCh3 <- 20 # 'south' (srtm) points
133
   
134
# NEW: Three plots: The difference between Mosaic and CDEM for pixels along
135
# each of three border edges: North (ASTER), Border, South (SRTM)
136
# for now, north, border, south have separate plots
137

    
138
# We need to tailor the plot 'triplet' X and Y axes to the dynamic ranges of all three data sets.
139

    
140
# Create values for the three 'delta across pixel boundaries' plots: 
141
#  1) Y-Axis: columnIDs 
142
#  2) Y-Axis: CDEM elevations
143
#  3) Y-Axis: 
144

    
145
   deltaBoundaryASTER <- northRowTbl$diffNorthMosaicCDEM - northRowTbl$diffSouthMosaicCDEM
146
   deltaBoundaryBorder <- borderRowTbl$diffNorthMosaicCDEM - borderRowTbl$diffSouthMosaicCDEM  
147
   deltaBoundarySRTM <- southRowTbl$diffNorthMosaicCDEM - southRowTbl$diffSouthMosaicCDEM   
148
   
149
   normDeltaBoundaryASTER <-  (deltaBoundaryASTER / northRowTbl$cdemNorth * 100)
150
   normDeltaBoundaryBorder <- (deltaBoundaryBorder / borderRowTbl$cdemNorth * 100)  
151
   normDeltaBoundarySRTM <-   (deltaBoundarySRTM / southRowTbl$cdemNorth * 100)
152
   
153
   par(mfrow=c(3,1)) # Create a set of three column multi-plots
154
   
155
   commonXAxis <- c(0,max(pointTable$ColumnID))   
156
   commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)   
157

    
158
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
159
   plot(northRowTbl$ColumnID,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
160
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
161
   plot(northRowTbl$ColumnID,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
162
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
163
   plot(northRowTbl$ColumnID,deltaBoundarySRTM,main="Boundary Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
164
message("ColumnID-X-Axis is done")
165
browser()
166
#dev.new()
167
#   par(mfrow=c(3,1)) # Create a three column multi-plot pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
168
   
169
   commonXAxis <- c(0,max(pointTable$cdemNorth))   
170
   commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)   
171
   
172
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
173
   plot(northRowTbl$cdemNorth,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
174
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
175
   plot(northRowTbl$cdemNorth,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
176
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
177
   plot(northRowTbl$cdemNorth,deltaBoundarySRTM,main="Boundary Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
178
message("ColumnID-CDEM Elevation done")
179
browser()
180
   commonXAxis <- c(0,max(pointTable$cdemNorth))   
181
   commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)   
182
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
183
   plot(northRowTbl$cdemNorth,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
184
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
185
   plot(northRowTbl$cdemNorth,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
186
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
187
   plot(northRowTbl$cdemNorth,normDeltaBoundarySRTM,main="Norm Bdry Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
188
message("NORMALIZED ColumnID-CDEM Elevation done")
189
browser()
190
   commonXAxis <- c(0,max(pointTable$ColumnID))   
191
   commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)   
192
   
193
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
194
   plot(northRowTbl$ColumnID,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
195
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
196
   plot(northRowTbl$ColumnID,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
197
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
198
   plot(northRowTbl$ColumnID,normDeltaBoundarySRTM,main="Norm Bdry Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
199
message("NORMALIZED ColumnID-X-Axis is done...")
200

    
201
# if 'plot flag' is set, send the plots to a JPG file. 
202

    
203
if (JpgPlotFileFlag)
204
{
205
  message("creating JPG plot files...")
206
  browser()
207
jpeg(file="RowBoundaryDeltaColIDXAxisASTER.jpg")
208
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
209
   plot(northRowTbl$ColumnID,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
210
dev.off()
211
jpeg(file="RowBoundaryDeltaColIDXAxisBORDER.jpg")   
212
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
213
   plot(northRowTbl$ColumnID,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
214
dev.off()
215
jpeg(file="RowBoundaryDeltaColIDXAxisSRTM.jpg")   
216
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
217
   plot(northRowTbl$ColumnID,deltaBoundarySRTM,main="Row Boundary Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
218
message("ColumnID-X-Axis plots are  done")
219
#dev.new()
220
   commonXAxis <- c(0,max(pointTable$cdemNorth))   
221
   commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)   
222
jpeg(file="RowBoundaryDeltaCDEMElevXAxisASTER.jpg")
223
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
224
   plot(northRowTbl$cdemNorth,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
225
dev.off()
226
jpeg(file="RowBoundaryDeltaCDEMElevXAxisBORDER.jpg")
227
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
228
   plot(northRowTbl$cdemNorth,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
229
dev.off()
230
jpeg(file="RowBoundaryDeltaCDEMElevXAxisCDEM.jpg")
231
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
232
   plot(northRowTbl$cdemNorth,deltaBoundarySRTM,main="Boundary Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
233
dev.off()
234
message("ColumnID-CDEM Elevation Plots are done")
235
#browser()
236
   commonXAxis <- c(0,max(pointTable$cdemNorth))   
237
   commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)   
238
jpeg(file="NormRowBoundaryDeltaColIDXAxisASTER.jpg")
239
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
240
   plot(northRowTbl$cdemNorth,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
241
dev.off()
242
jpeg(file="NormRowBoundaryDeltaColIDXAxisBORDER.jpg")
243
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
244
   plot(northRowTbl$cdemNorth,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
245
dev.off()
246
jpeg(file="NormRowBoundaryDeltaColIDXAxisSRTM.jpg")
247
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
248
   plot(northRowTbl$cdemNorth,normDeltaBoundarySRTM,main="Norm Bdry Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
249
dev.off()
250
message("NORMALIZED ColumnID-CDEM Elevation Plots are done")
251
browser()
252
   commonXAxis <- c(0,max(pointTable$ColumnID))   
253
   commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)   
254
jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisASTER.jpg")   
255
   xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
256
   plot(northRowTbl$ColumnID,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
257
dev.off()
258
jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisBORDER.jpg")
259
   xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
260
   plot(northRowTbl$ColumnID,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
261
dev.off()
262
jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisSRTM.jpg")
263
   xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
264
   plot(northRowTbl$ColumnID,normDeltaBoundarySRTM,main="Norm Bdry Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
265
dev.off()   
266
message("NORMALIZED ColumnID-X-Axis Plot files have been created..")
267
}
268
message("...All plots created - hit key to delete them...")
269
browser()
270
graphics.off()
271
} 
(2-2/13)