##############################################################################################

#

# CreateImageDiffPlots

# 

# R script generates collection plots that explore elevation discontinuities 

# border between ASTER and STRM data in mosaic images under development for Iplant.

# Plots overwritten to single graphics device for first version. 

# Panel Plots to be added to next version.

#

# To run:

#

#  2) start R

#  3) > source("MosaicCDEMDifferencePlots.r")

# Input file: Csv table containing elevations for randomly-selected point pairs from 

# ASTER/CGIAR mosaic and CDEM images. Table is generated by the script: 'makeImagePairTable.r'

# Author: Rick Reeves, NCEAS

# May 1, 2011

##############################################################################################

#CreateImageDiffPlots <- function()

#{

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

#

# Table created by randomly sampling from two superimposed 

# image: 

#  1) DOM mosaic image comprised of ASTER and SRTM components.

#  2) 'Baseline' Canadian DEM (CDEM) image.

# 

# The input table contains three rows for each randomly-selected

# pixel pair from both images. Each row contains two pixel pairs,

# the first pair drawn from the image mosaic, the second pair

# drawn from the CDEM image: 

#   First pair: North pixel, South pixel (ASTER/SRTM mosaic)

#   Second pair: North pixel, South pixel (CDEM)

#

#  The first row of each 'triplet' contains pixel pairs North of border,

#  The second row contains pixel pairs spanning  border,

#  The third row contains pixel pairs South of border,

#

# This script generates a series of plots that display

# differences between:

#    1) The mosaic and CDEM images

#    2) Image pixels on and away from the ASTER / SRTM boundary.

# 

   northRowIndexes = seq(from=1, to=(nrow(pointTable) - 3),by=3)

   borderRowIndexes = seq(from=2, to=(nrow(pointTable) - 1),by=3)   

   southRowIndexes = seq(from=3, to=(nrow(pointTable)),by=3)

#

# calculate and append the difference between elevations

# and CDEM

# these lines create the inputs for differnce image plots for each of three

# pixel pair subsets: North of border (All Aster), border (combo Aster/Srtm),

#                     South of border (all Srtm)

# 

# First, add the 'difference columns' to the entire table

#

   pointTable <-cbind(pointTable,(pointTable$elevNorth - pointTable$elevSouth))

   pointTable <-cbind(pointTable,(pointTable$cdemNorth - pointTable$cdemSouth))

   pointTable <-cbind(pointTable,(pointTable$elevNorth - pointTable$cdemNorth))

   pointTable <-cbind(pointTable,(pointTable$elevSouth - pointTable$cdemSouth))

#   

   colnames(pointTable)[6] <- "diffMosaicNorthSouth"

   colnames(pointTable)[7] <- "diffCDEMNorthSouth"

   colnames(pointTable)[8] <- "diffNorthMosaicCDEM"

   colnames(pointTable)[9] <- "diffSouthMosaicCDEM"   

# Difference between Mosaic (ASTER or CGIAR or border) 

# and CDEM elevation as pertentage of the mosaic elevation

   pointTable <-cbind(pointTable,(pointTable$diffSouthMosaicCDEM/pointTable$elevSouth * 100)) 

   colnames(pointTable)[10] <- "magDiffMosaicCDEMNorthPct"

   colnames(pointTable)[11] <- "magDiffMosaicCDEMSouthPct"   

# For the plots, subdivide the table into three segments: 

#     rows north of border

#     rows crossing border

#     rows south of border

   northRowTbl = pointTable[northRowIndexes,]

   borderRowTbl = pointTable[borderRowIndexes,]

   southRowTbl = pointTable[southRowIndexes,]

# Create plots

#

# First, plot North/South pixel elevation for North/Border/South subsets 'distance East from Western edge'.

#

   plot(northRowTbl$ColumnID,northRowTbl$elevNorth,main="North (ASTER) Mosaic North Pixel Elevation",col="red",pch=".")

browser()   

browser()   

browser()   

browser()   

browser()   

browser()

browser()   

browser()   

browser()

browser()   

browser()   

browser()   

# Next, (for now) difference between North and South pixels in pair, for various regions

# 

   plot(northRowTbl$ColumnID,northRowTbl$diffMosaicNorthSouth,main="North (ASTER) Mosaic North/South Diff",col="red",pch=".")

browser()   

browser()   

browser()   

   plot(northRowTbl$ColumnID,northRowTbl$diffCDEMNorthSouth,main="North (ASTER) CDEM North/South Diff",col="red",pch=".")

browser()   

browser()   

browser()   

# Next: the Mosaic/CDEM difference divided by the (ASTER or SRTM) elevation

browser()  

browser()  

browser()  

  plot(northRowTbl$elevNorth,northRowTbl$diffCDEMNorthSouth,col="red",pch=".",main="ASTER Elev Diff vs Mosaic Elev (North Pixel)")

# plot(northRowTbl$elevSouth,southRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="ASTER Elev Diff vs Mosaic Elev (South Pixel)")  

browser()  

#  plot(borderRowTbl$elevSouth,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="Border Elev Diff vs Mosaic Elev (South Pixel)")  

browser()  

#  plot(southRowTbl$elevSouth,southRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="SRTM Elev Diff vs Mosaic Elev (South Pixel)")  

#}