/ - Diff - Environment and organisms - NCEAS Projects

« Previous | Next »

Revision a09a09c0

Added by Jim Regetz over 12 years ago

ID a09a09c0b9afeb1e0aa655e4ebd758eb5cdcfdd5

removed obsolete/preliminary DEM 60N assessment scripts

     ##############################################################################################
+    #
     # CreateImageDiffPlots
+    #
     # R script generates collection of nine plots that display the distributions of  populations of
     # three elevation pixel pairs in proximity to the border between ASTER and STRM data in mosaic
     # images# under development for the Environment and Organisms project.
+    #
+    #
     # Inputs:
     #   1) Comma Separated Value (CSV) table containing randomly-sampled elevation value pairs,
     #      extracted from ASTER/CGIAR mosaic and CDEM images, created by the R script: makeImagePairTable.r
     #      Note: At present, be sure that this file has been sorted by column 1 (ColumnID) in Excel
     #      before using in this program.
+    #
     #   2) sampling factor: integer (range=1 to number of recors in input table. Determines
     #      the size of randomly-selected sampe of total table record 'triplets' (north, border,
     #      and south) rows of each column sample) to be displayed in the plots:
     #           sampling factor = 1  : plot every table record
     #                             10 : plot a random sample containing 1/10th of records
     #                            100 : plot random sample containing 1/100th of records.
+    #
     # To run:
+    #
     #  1) place this file and the input file "tableForMark4000_5_8_SortColID.csv" in folder
     #  2) start R
     #  3) > source("CreateImageDiffPlotsOverlay.r")
     #  4) > CreateImageDiffPlots(sampling factor)
     #     < press <cr> key to view plots sequentially
+    #
     # TODO: add symbol legend to each plot.
     # Author: Rick Reeves, NCEAS
     # May 14, 2011
     ##############################################################################################
     CreateImageDiffPlots <- function(plotSampFact = 100)
+    {
     # Check plotSampleFact range
        if (plotSampFact < 1)
           plotSampFact = 1
        if (plotSampFact > 100)
           plotSampFact = 100
     # Read input table that was sorted OFFLINE in Excel on the first column (ColumnID)
        pointTable <-read.csv("tableForMark4000_5_8_SortColID.csv")
     # 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$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
        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
        northRowTblAll = pointTable[northRowIndexes,]
        borderRowTblAll = pointTable[borderRowIndexes,]
        southRowTblAll = pointTable[southRowIndexes,]
        subset <- 1:nrow(northRowTblAll)
        randSub <- sample(subset,as.integer(length(subset) / plotSampFact))
        northRowTbl <- northRowTblAll[randSub,]
        borderRowTbl <- borderRowTblAll[randSub,]
        southRowTbl <- southRowTblAll[randSub,]
     message("hit key to create each plot...")
     browser()
     # Three plotting characters used
        plotCh1 <- 17 # 'north' (Aster) points
        plotCh2 <- 18 # 'border' (Aster+srtm) points
        plotCh3 <- 20 # 'south' (srtm) points
     # First plot pair: North/South Mosaic Image pixel elevation for North/South subsets 'distance East from Western edge'.
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$elevNorth,na.rm=TRUE),sd(northRowTbl$elevNorth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$elevNorth,na.rm=TRUE),sd(borderRowTbl$elevNorth,na.rm=TRUE),
                                                                                 mean(southRowTbl$elevNorth,na.rm=TRUE),sd(southRowTbl$elevNorth,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$elevNorth,main="North Mosaic Pixel Elev: ASTER (red) Border (green) SRTM (blue)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$elevNorth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$elevNorth,col="blue",pch=plotCh3)
     browser()
     dev.new()
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$elevSouth,na.rm=TRUE),sd(northRowTbl$elevSouth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$elevSouth,na.rm=TRUE),sd(borderRowTbl$elevSouth,na.rm=TRUE),
                                                                                 mean(southRowTbl$elevSouth,na.rm=TRUE),sd(southRowTbl$elevSouth,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$elevSouth,main="South Mosaic Pixel Elev: ASTER (red) Border (green) SRTM (blue)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$elevSouth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$elevSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
     # Second plot pair: North/South CDEM pixel elevation for North/South subsets 'distance East from Western edge'.
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$cdemNorth,na.rm=TRUE),sd(northRowTbl$cdemNorth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$cdemNorth,na.rm=TRUE),sd(borderRowTbl$cdemNorth,na.rm=TRUE),
                                                                                 mean(southRowTbl$cdemNorth,na.rm=TRUE),sd(southRowTbl$cdemNorth,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$cdemNorth,main="North CDEM Pixel Elev ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$cdemNorth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$cdemNorth,col="blue",pch=plotCh3)
     browser()
     dev.new()
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$cdemSouth,na.rm=TRUE),sd(northRowTbl$cdemSouth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$cdemSouth,na.rm=TRUE),sd(borderRowTbl$cdemSouth,na.rm=TRUE),
                                                                                 mean(southRowTbl$cdemSouth,na.rm=TRUE),sd(southRowTbl$cdemSouth,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$cdemSouth,main="South CDEM Pixel Elev ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$cdemSouth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$cdemSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
     # Third plot pair: difference between North and South pixels in pair.
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$diffMosaicNorthSouth,na.rm=TRUE),sd(northRowTbl$diffMosaicNorthSouth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$diffMosaicNorthSouth,na.rm=TRUE),sd(borderRowTbl$diffMosaicNorthSouth,na.rm=TRUE),
                                                                                 mean(southRowTbl$diffMosaicNorthSouth,na.rm=TRUE),sd(southRowTbl$diffMosaicNorthSouth,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$diffMosaicNorthSouth,main="Mosaic N/S Pixel Diff: ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$diffMosaicNorthSouth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$diffMosaicNorthSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$diffCDEMNorthSouth,na.rm=TRUE),sd(northRowTbl$diffCDEMNorthSouth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$diffCDEMNorthSouth,na.rm=TRUE),sd(borderRowTbl$diffCDEMNorthSouth,na.rm=TRUE),
                                                                                 mean(southRowTbl$diffCDEMNorthSouth,na.rm=TRUE),sd(southRowTbl$diffCDEMNorthSouth,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$diffCDEMNorthSouth,main="CDEM N/S Pixel Diff: ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$diffCDEMNorthSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
     # Fourth plot pair: the Mosaic/CDEM difference for North and South pixels
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$diffNorthMosaicCDEM,na.rm=TRUE),sd(northRowTbl$diffNorthMosaicCDEM,na.rm=TRUE),
                                                                                 mean(borderRowTbl$diffNorthMosaicCDEM,na.rm=TRUE),sd(borderRowTbl$diffNorthMosaicCDEM,na.rm=TRUE),
                                                                                 mean(southRowTbl$diffNorthMosaicCDEM,na.rm=TRUE),sd(southRowTbl$diffNorthMosaicCDEM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$diffMosaicNorthSouth,main="Mosaic / CDEM Pixel Diff (North): ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$diffMosaicNorthSouth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$diffMosaicNorthSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
        xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$diffSouthMosaicCDEM,na.rm=TRUE),sd(northRowTbl$diffSouthMosaicCDEM,na.rm=TRUE),
                                                                                 mean(borderRowTbl$diffSouthMosaicCDEM,na.rm=TRUE),sd(borderRowTbl$diffSouthMosaicCDEM,na.rm=TRUE),
                                                                                 mean(southRowTbl$diffSouthMosaicCDEM,na.rm=TRUE),sd(southRowTbl$diffSouthMosaicCDEM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,northRowTbl$diffCDEMNorthSouth,main="Mosaic / CDEM Pixel Diff (South): ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
        points(borderRowTbl$ColumnID,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=plotCh2)
        points(southRowTbl$ColumnID,southRowTbl$diffCDEMNorthSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
     # Fifth plot pair: the Mosaic/CDEM difference divided by the (ASTER or SRTM) elevation
       xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),sd(northRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),
                                                                                 mean(borderRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),sd(borderRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),
                                                                                 mean(southRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),sd(southRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE))
       plot(northRowTbl$ColumnID,northRowTbl$magDiffMosaicCDEMNorthPct,main="Pixel Elev Diff as % of Mosaic Elev (North) - ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
       points(borderRowTbl$ColumnID,borderRowTbl$magDiffMosaicCDEMNorthPct,col="darkgreen",pch=plotCh2)
       points(southRowTbl$ColumnID,southRowTbl$magDiffMosaicCDEMNorthPct,col="blue",pch=plotCh3)
     browser()
     dev.new()
       xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$magDiffMosaicCDEMSouthPct,na.rm=TRUE),sd(northRowTbl$magDiffMosaicCDEMSouthPct,na.rm=TRUE),
                                                                                 mean(borderRowTbl$magDiffMosaicCDEMSouthPct,na.rm=TRUE),sd(borderRowTbl$magDiffMosaicCDEMSouthPct,na.rm=TRUE),
                                                                                  mean(southRowTbl$magDiffMosaicCDEMSouthPct,na.rm=TRUE),sd(southRowTbl$magDiffMosaicCDEMSouthPct,na.rm=TRUE))
       plot(northRowTbl$ColumnID,northRowTbl$magDiffMosaicCDEMSouthPct,main="Pixel Elev Diff as % of Mosaic Elev (South) - ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
       points(borderRowTbl$ColumnID,borderRowTbl$magDiffMosaicCDEMSouthPct,col="darkgreen",pch=plotCh2)
       points(southRowTbl$ColumnID,southRowTbl$magDiffMosaicCDEMSouthPct,col="blue",pch=plotCh3)
     browser()
     dev.new()
     # Final plot pair: Mosaic vs CDEM difference vs elevation (north or south)
     xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$diffMosaicNorthSouth,na.rm=TRUE),sd(northRowTbl$diffMosaicNorthSouth,na.rm=TRUE),
                                                                               mean(borderRowTbl$diffMosaicNorthSouth,na.rm=TRUE),sd(borderRowTbl$diffMosaicNorthSouth,na.rm=TRUE),
                                                                               mean(southRowTbl$diffMosaicNorthSouth,na.rm=TRUE),sd(southRowTbl$diffMosaicNorthSouth,na.rm=TRUE))
     plot(northRowTbl$elevNorth,northRowTbl$diffMosaicNorthSouth,main="Mosaic Elev N/S Diff ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
     points(borderRowTbl$elevNorth,borderRowTbl$diffMosaicNorthSouth,col="darkgreen",pch=plotCh2)
     points(southRowTbl$elevNorth,southRowTbl$diffMosaicNorthSouth,col="blue",pch=plotCh3)
     browser()
     dev.new()
       xAxisLbl <- sprintf("M/SD: AST: %.2f / %.2f BRD: %.2f / %.2f STM: %.2f / %.2f",mean(northRowTbl$diffCDEMNorthSouth,na.rm=TRUE),sd(northRowTbl$diffCDEMNorthSouth,na.rm=TRUE),
                                                                                 mean(borderRowTbl$diffCDEMNorthSouth,na.rm=TRUE),sd(borderRowTbl$diffCDEMNorthSouth,na.rm=TRUE),
                                                                                 mean(southRowTbl$diffCDEMNorthSouth,na.rm=TRUE),sd(southRowTbl$diffCDEMNorthSouth,na.rm=TRUE))
       plot(northRowTbl$elevNorth,northRowTbl$diffCDEMNorthSouth,main="CDEM Elev N/S Diff ASTER (r) Border (g) SRTM (b)",xlab=xAxisLbl,col="red",pch=plotCh1)
       points(borderRowTbl$elevNorth,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=plotCh2)
       points(southRowTbl$elevNorth,southRowTbl$diffCDEMNorthSouth,col="blue",pch=plotCh3)
     message("All plots created - hit key to delete them...")
     browser()
     graphics.off()
+    }

     ##############################################################################################
+    #
     # CreateImageDiffPlots
+    #
     # R script generates collection of nine plots that display the distributions of  populations of
     # three elevation pixel pairs in proximity to the border between ASTER and STRM data in mosaic
     # images# under development for the Environment and Organisms project.
+    #
+    #
     # Inputs:
     #   1) Comma Separated Value (CSV) table containing randomly-sampled elevation value pairs,
     #      extracted from ASTER/CGIAR mosaic and CDEM images, created by the R script: makeImagePairTable.r
     #      Note: At present, be sure that this file has been sorted by column 1 (ColumnID) in Excel
     #      before using in this program.
     #      NOTE: This filename is set in the body of the script.
+    #
     #   2) sampling factor: integer (range=1 to number of recors in input table. Determines
     #      the size of randomly-selected sampe of total table record 'triplets' (north, border,
     #      and south) rows of each column sample) to be displayed in the plots:
     #           sampling factor = 1  : plot every table record
     #                             10 : plot a random sample containing 1/10th of records
     #                            100 : plot random sample containing 1/100th of records.
     #   3) JpgPlotFileFlag: Set to TRUE to have *individual* plots written to JPG files.
+    #
     # To run:
+    #
     #  1) place this file and the input file "tableForMark4000_5_8_SortColID.csv" in folder
     #  2) start R
     #  3) > source("CreateImageDiffPlotsOverlayNewPlots.r")
     #  4) > CreateImageDiffPlots(sampling factor, JpegFlag)
     #     < press <cr> key to view plots sequentially
+    #
     # TODO: add symbol legend to each plot.
     # Author: Rick Reeves, NCEAS
     # May 14, 2011
     # May 17, 2011: This version generates 'delta scatterplots' specified by Mark and Jim.
     # May 18, 2011: If JpgPlotFlag set, write individual plot files to JPG files.
     ##############################################################################################
     CreateImageDiffPlots <- function(plotSampFact = 1,JpgPlotFileFlag = TRUE)
+    {
     # Check plotSampleFact range
        if (plotSampFact < 1)
           plotSampFact = 1
     # Read input table that was sorted OFFLINE in Excel on the first column (ColumnID)
        pointTable <-read.csv("pixelPairs36000_5_8_TS.csv")
     # 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"
     # add a placeholder for the 'boundary' value, across each table
     #  pointTable <-cbind(pointTable,(-1))
     #  colnames(pointTable)[10] <- "deltaAcrossBorder"
     # Difference between Mosaic (ASTER or CGIAR or border)
     # and CDEM elevation as pertentage of the mosaic elevation
        pointTable <-cbind(pointTable,(pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
        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
        northRowTblAll = pointTable[northRowIndexes,]
        borderRowTblAll = pointTable[borderRowIndexes,]
        southRowTblAll = pointTable[southRowIndexes,]
        subset <- 1:nrow(northRowTblAll)
        randSub <- sample(subset,as.integer(length(subset) / plotSampFact))
        northRowTbl <- northRowTblAll[randSub,]
        borderRowTbl <- borderRowTblAll[randSub,]
        southRowTbl <- southRowTblAll[randSub,]
     message("hit key to create each plot...")
     browser()
     # Three plotting characters used
        plotCh1 <- 17 # 'north' (aster) points
        plotCh2 <- 18 # 'border' (aster+srtm) points
        plotCh3 <- 20 # 'south' (srtm) points
     # NEW: Three plots: The difference between Mosaic and CDEM for pixels along
     # each of three border edges: North (ASTER), Border, South (SRTM)
     # for now, north, border, south have separate plots
     # We need to tailor the plot 'triplet' X and Y axes to the dynamic ranges of all three data sets.
     # Create values for the three 'delta across pixel boundaries' plots:
     #  1) Y-Axis: columnIDs
     #  2) Y-Axis: CDEM elevations
     #  3) Y-Axis:
        deltaBoundaryASTER <- northRowTbl$diffNorthMosaicCDEM - northRowTbl$diffSouthMosaicCDEM
        deltaBoundaryBorder <- borderRowTbl$diffNorthMosaicCDEM - borderRowTbl$diffSouthMosaicCDEM
        deltaBoundarySRTM <- southRowTbl$diffNorthMosaicCDEM - southRowTbl$diffSouthMosaicCDEM
        normDeltaBoundaryASTER <-  (deltaBoundaryASTER / northRowTbl$cdemNorth * 100)
        normDeltaBoundaryBorder <- (deltaBoundaryBorder / borderRowTbl$cdemNorth * 100)
        normDeltaBoundarySRTM <-   (deltaBoundarySRTM / southRowTbl$cdemNorth * 100)
        par(mfrow=c(3,1)) # Create a set of three column multi-plots
        commonXAxis <- c(0,max(pointTable$ColumnID))
        commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundarySRTM,main="Boundary Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     message("ColumnID-X-Axis is done")
     browser()
     #dev.new()
     #   par(mfrow=c(3,1)) # Create a three column multi-plot pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
        commonXAxis <- c(0,max(pointTable$cdemNorth))
        commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundarySRTM,main="Boundary Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     message("ColumnID-CDEM Elevation done")
     browser()
        commonXAxis <- c(0,max(pointTable$cdemNorth))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundarySRTM,main="Norm Bdry Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     message("NORMALIZED ColumnID-CDEM Elevation done")
     browser()
        commonXAxis <- c(0,max(pointTable$ColumnID))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundarySRTM,main="Norm Bdry Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     message("NORMALIZED ColumnID-X-Axis is done...")
     # if 'plot flag' is set, send the plots to a JPG file.
     if (JpgPlotFileFlag)
+    {
       message("creating JPG plot files...")
       browser()
     jpeg(file="RowBoundaryDeltaColIDXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="RowBoundaryDeltaColIDXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="RowBoundaryDeltaColIDXAxisSRTM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundarySRTM,main="Row Boundary Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     message("ColumnID-X-Axis plots are  done")
     #dev.new()
        commonXAxis <- c(0,max(pointTable$cdemNorth))
        commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)
     jpeg(file="RowBoundaryDeltaCDEMElevXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="RowBoundaryDeltaCDEMElevXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="RowBoundaryDeltaCDEMElevXAxisCDEM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundarySRTM,main="Boundary Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     dev.off()
     message("ColumnID-CDEM Elevation Plots are done")
     #browser()
        commonXAxis <- c(0,max(pointTable$cdemNorth))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
     jpeg(file="NormRowBoundaryDeltaColIDXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaColIDXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaColIDXAxisSRTM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundarySRTM,main="Norm Bdry Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     dev.off()
     message("NORMALIZED ColumnID-CDEM Elevation Plots are done")
     browser()
        commonXAxis <- c(0,max(pointTable$ColumnID))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
     jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisSRTM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundarySRTM,main="Norm Bdry Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     dev.off()
     message("NORMALIZED ColumnID-X-Axis Plot files have been created..")
+    }
     message("...All plots created - hit key to delete them...")
     browser()
     graphics.off()
+    }

     ##############################################################################################
+    #
     # CreateImageDiffPlots
+    #
     # R script generates collection of nine plots that display the distributions of  populations of
     # three elevation pixel pairs in proximity to the border between ASTER and STRM data in mosaic
     # images# under development for the Environment and Organisms project.
+    #
+    #
     # Inputs:
     #   1) Comma Separated Value (CSV) table containing randomly-sampled elevation value pairs,
     #      extracted from ASTER/CGIAR mosaic and CDEM images, created by the R script: makeImagePairTable.r
     #      Note: At present, be sure that this file has been sorted by column 1 (ColumnID) in Excel
     #      before using in this program.
     #      NOTE: This filename is set in the body of the script.
+    #
     #   2) sampling factor: integer (range=1 to number of recors in input table. Determines
     #      the size of randomly-selected sampe of total table record 'triplets' (north, border,
     #      and south) rows of each column sample) to be displayed in the plots:
     #           sampling factor = 1  : plot every table record
     #                             10 : plot a random sample containing 1/10th of records
     #                            100 : plot random sample containing 1/100th of records.
     #   3) JpgPlotFileFlag: Set to TRUE to have *individual* plots written to JPG files.
+    #
     # To run:
+    #
     #  1) place this file and the input file "tableForMark4000_5_8_SortColID.csv" in folder
     #  2) start R
     #  3) > source("CreateImageDiffPlotsOverlay.r")
     #  4) > CreateImageDiffPlots(sampling factor,JpegFlag)
     #     < press <cr> key to view plots sequentially
+    #
     # TODO: add symbol legend to each plot.
     # Author: Rick Reeves, NCEAS
     # May 14, 2011
     # May 17, 2011: This version generates 'delta scatterplots' specified by Mark and Jim.
     # May 18, 2011: If JpgPlotFlag set, write individual plot files to JPG files.
     # May 22, 2011: add Lowess lines to center of scatterplot distribution
     ##############################################################################################
     CreateImageDiffPlotsNewLL <- function(plotSampFact = 1,JpgPlotFileFlag = TRUE)
+    {
     # Check plotSampleFact rangels -l *.tif
        if (plotSampFact < 1)
           plotSampFact = 1
     # Read input table that was sorted OFFLINE in Excel on the first column (ColumnID)
        pointTable <-read.csv("pixelPairs36000_5_8EvenSortCol1.csv")
     #   pointTable <-read.csv("pixelPairs36000_Exa.csv")
     # 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"
     # add a placeholder for the 'boundary' value, across each table
     #  pointTable <-cbind(pointTable,(-1))
     #  colnames(pointTable)[10] <- "deltaAcrossBorder"
     # Difference between Mosaic (ASTER or CGIAR or border)
     # and CDEM elevation as pertentage of the mosaic elevation
        pointTable <-cbind(pointTable,(pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
        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
        northRowTblAll = pointTable[northRowIndexes,]
        borderRowTblAll = pointTable[borderRowIndexes,]
        southRowTblAll = pointTable[southRowIndexes,]
        subset <- 1:nrow(northRowTblAll)
        randSub <- sample(subset,as.integer(length(subset) / plotSampFact))
        northRowTbl <- northRowTblAll[randSub,]
        borderRowTbl <- borderRowTblAll[randSub,]
        southRowTbl <- southRowTblAll[randSub,]
     message("hit key to create each plot...")
     browser()
     # Three plotting characters used
        plotCh1 <- 17 # 'north' (aster) points
        plotCh2 <- 18 # 'border' (aster+srtm) points
        plotCh3 <- 20 # 'south' (srtm) points
        lowessLineColor <- "yellow"
     # NEW: Three plots: The difference between Mosaic and CDEM for pixels along
     # each of three border edges: North (ASTER), Border, South (SRTM)
     # for now, north, border, south have separate plots
     # We need to tailor the plot 'triplet' X and Y axes to the dynamic ranges of all three data sets.
     # Create values for the three 'delta across pixel boundaries' plots:
     #  1) Y-Axis: columnIDs
     #  2) Y-Axis: CDEM elevations
     #  3) Y-Axis:
        deltaBoundaryASTER <- northRowTbl$diffNorthMosaicCDEM - northRowTbl$diffSouthMosaicCDEM
        deltaBoundaryBorder <- borderRowTbl$diffNorthMosaicCDEM - borderRowTbl$diffSouthMosaicCDEM
        deltaBoundarySRTM <- southRowTbl$diffNorthMosaicCDEM - southRowTbl$diffSouthMosaicCDEM
        normDeltaBoundaryASTER <-  (deltaBoundaryASTER / northRowTbl$cdemNorth * 100)
        normDeltaBoundaryBorder <- (deltaBoundaryBorder / borderRowTbl$cdemNorth * 100)
        normDeltaBoundarySRTM <-   (deltaBoundarySRTM / southRowTbl$cdemNorth * 100)
        par(mfrow=c(1,3)) # Create a set of three column multi-plots
        commonXAxis <- c(0,max(pointTable$ColumnID))
        commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Delta: ASTER - CDEM",xlab=xAxisLbl,col="red",pch=plotCh1)
        vv <- cbind(northRowTbl$ColumnID,deltaBoundaryASTER)
     #   vx <- vv[-(which(is.na(vv[,"deltaBoundaryASTER"]))),]
        vx <- vv[(which(!is.na(vv[,"deltaBoundaryASTER"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Delta: Bdry Mix - CDEM",sub="X Axis: West-East Columns",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        vv <- cbind(northRowTbl$ColumnID,deltaBoundaryBorder)
        vx <- vv[(which(!is.na(vv[,"deltaBoundaryBorder"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundarySRTM,main="Delta: SRTM - CDEM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
        vv <- cbind(northRowTbl$ColumnID,deltaBoundarySRTM)
        vx <- vv[(which(!is.na(vv[,"deltaBoundarySRTM"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
     message("ColumnID-X-Axis is done")
     browser()
     #dev.new()
     #   par(mfrow=c(3,1)) # Create a three column multi-plot pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
        commonXAxis <- c(0,max(pointTable$cdemNorth,na.rm=TRUE))
        commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)
     #   northDelta.lo <- loess(deltaBoundaryASTER ~ deltaBoundaryASTER,northRowTbl)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Delta: ASTER - CDEM",xlab=xAxisLbl,col="red",pch=plotCh1)
        vv <- cbind(northRowTbl$cdemNorth,deltaBoundaryASTER)
        vx <- vv[(which(!is.na(vv[,"deltaBoundaryASTER"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main=" Delta: Bdy Mix - CDEM",sub="X Axis: CDEM Elevation",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        vv <- cbind(northRowTbl$cdemNorth,deltaBoundaryBorder)
        vx <- vv[(which(!is.na(vv[,"deltaBoundaryBorder"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundarySRTM,main="Delta: SRTM - CDEM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
        vv <- cbind(northRowTbl$cdemNorth,deltaBoundarySRTM)
        vx <- vv[(which(!is.na(vv[,"deltaBoundarySRTM"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
     message("ColumnID-CDEM Elevation done")
     browser()
        commonXAxis <- c(0,max(pointTable$cdemNorth,na.rm=TRUE))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Delta Norm Elev: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
        vv <- cbind(northRowTbl$cdemNorth,normDeltaBoundaryASTER)
        vx <- vv[(which(!is.na(vv[,"normDeltaBoundaryASTER"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Delta Norm Elev: Border",sub="X Axis: CDEM Elevation",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        vv <- cbind(northRowTbl$cdemNorth,normDeltaBoundaryBorder)
        vx <- vv[(which(!is.na(vv[,"normDeltaBoundaryBorder"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundarySRTM,main="Delta Norm Elev: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
        vv <- cbind(northRowTbl$cdemNorth,normDeltaBoundarySRTM)
        vx <- vv[(which(!is.na(vv[,"normDeltaBoundarySRTM"]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
     message("NORMALIZED ColumnID-CDEM Elevation done")
     browser()
        commonXAxis <- c(0,max(pointTable$cdemNorth,na.rm=TRUE))
        commonYAxis <- range(c(northRowTbl$magDiffMosaicCDEMNorthPct,borderRowTbl$magDiffMosaicCDEMNorthPct,southRowTbl$magDiffMosaicCDEMNorthPct),na.rm=TRUE)
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(northRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),sd(northRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,northRowTbl$magDiffMosaicCDEMNorthPct,xlim=commonXAxis, ylim=commonYAxis,main="Delta Vs Elev Mag: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
        vv <- cbind(northRowTbl$cdemNorth,northRowTbl$magDiffMosaicCDEMNorthPct)
        vx <- vv[(which(!is.na(vv[,2]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(borderRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),sd(borderRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,borderRowTbl$magDiffMosaicCDEMNorthPct,xlim=commonXAxis, ylim=commonYAxis,main="Delta Vs Elev Mag:: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
        vv <- cbind(northRowTbl$cdemNorth,deltaBoundarySRTM)
        vx <- vv[(which(!is.na(vv[,2]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(southRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE),sd(southRowTbl$magDiffMosaicCDEMNorthPct,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,southRowTbl$magDiffMosaicCDEMNorthPct,main="Delta Vs Elev Mag:: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
        vv <- cbind(northRowTbl$cdemNorth,southRowTbl$magDiffMosaicCDEMNorthPct)
        vx <- vv[(which(!is.na(vv[,2]))),]
        lines(stats::lowess(vx),col=lowessLineColor)
     message("Delta as Fcn Elev Mag-X-Axis is done...")
     # if 'plot flag' is set, send the plots to a JPG file.
     if (JpgPlotFileFlag)
+    {
       message("creating JPG plot files...")
       browser()
     jpeg(file="RowBoundaryDeltaColIDXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="RowBoundaryDeltaColIDXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Row Boundary Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="RowBoundaryDeltaColIDXAxisSRTM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,deltaBoundarySRTM,main="Row Boundary Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     message("ColumnID-X-Axis plots are  done")
     #dev.new()
        commonXAxis <- c(0,max(pointTable$cdemNorth))
        commonYAxis <- range(c(deltaBoundarySRTM,deltaBoundaryBorder,deltaBoundaryASTER),na.rm=TRUE)
     jpeg(file="RowBoundaryDeltaCDEMElevXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(deltaBoundaryASTER,na.rm=TRUE),sd(deltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="RowBoundaryDeltaCDEMElevXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(deltaBoundaryBorder,na.rm=TRUE),sd(deltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Boundary Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="RowBoundaryDeltaCDEMElevXAxisCDEM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(deltaBoundarySRTM,na.rm=TRUE),sd(deltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,deltaBoundarySRTM,main="Boundary Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     dev.off()
     message("ColumnID-CDEM Elevation Plots are done")
     #browser()
        commonXAxis <- c(0,max(pointTable$cdemNorth))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
     jpeg(file="NormRowBoundaryDeltaColIDXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaColIDXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta (CDEM Elev): Border",sub="vs Elev",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaColIDXAxisSRTM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$cdemNorth,normDeltaBoundarySRTM,main="Norm Bdry Delta (CDEM Elev): SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     dev.off()
     message("NORMALIZED ColumnID-CDEM Elevation Plots are done")
     browser()
        commonXAxis <- c(0,max(pointTable$ColumnID))
        commonYAxis <- range(c(normDeltaBoundarySRTM,normDeltaBoundaryBorder,normDeltaBoundaryASTER),na.rm=TRUE)
     jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisASTER.jpg")
        xAxisLbl <- sprintf("M/SD: ASTER: %.2f / %.2f",mean(normDeltaBoundaryASTER,na.rm=TRUE),sd(normDeltaBoundaryASTER,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundaryASTER,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: ASTER",xlab=xAxisLbl,col="red",pch=plotCh1)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisBORDER.jpg")
        xAxisLbl <- sprintf("M/SD: BORDER: %.2f / %.2f",mean(normDeltaBoundaryBorder,na.rm=TRUE),sd(normDeltaBoundaryBorder,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundaryBorder,xlim=commonXAxis, ylim=commonYAxis,main="Norm Bdry Delta: Border",sub="E-W Col",xlab=xAxisLbl,col="darkgreen",pch=plotCh2)
     dev.off()
     jpeg(file="NormRowBoundaryDeltaCDEMElevXAxisSRTM.jpg")
        xAxisLbl <- sprintf("M/SD: SRTM: %.2f / %.2f",mean(normDeltaBoundarySRTM,na.rm=TRUE),sd(normDeltaBoundarySRTM,na.rm=TRUE))
        plot(northRowTbl$ColumnID,normDeltaBoundarySRTM,main="Norm Bdry Delta: SRTM",xlim=commonXAxis, ylim=commonYAxis,xlab=xAxisLbl,col="blue",pch=plotCh3)
     dev.off()
     message("NORMALIZED ColumnID-X-Axis Plot files have been created..")
+    }
     message("...All plots created - hit key to delete them...")
     browser()
     graphics.off()
+    }

     # DEMBdyFixExpDecay - Implements Exponential Decay method to 'blend' DEM boundary edges
+    #
     # General form: decayFactor = e(-(d * DKc))
     #                 where d = distance interval from origin (in this case, an image row index)
     #                       DKc = empirically derived constant
+    #
+    #
     # Inputs:
     #   1) Blended Aster/SRTM DEM mosaic image.
     #   2) 'SRTM-Only' subimage corresponding to the NORTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the SRTM component of the mosaic.
     #   3) 'Aster-Only' subimage corresponding to the SOUTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the ASTER component of the mosaic.
     #   4) Other parameters, as they are needed
+    #
     # To run:
     #  2) > source("DEMBdyFixExpDecay.r")
     #  3) > DEMBdyFixExpDecay()
+    #
     # Author: Rick Reeves, NCEAS
     # June, 2011
     ##############################################################################################
     DEMBdyFixExpDecay <- function()
+    {
       require(raster)
       require(rgdal)
     # Read the mosaic image and component sub-images:
     #   - Aster/SRTM mosaic image (with boundary edge to be fixed)
     #   - 'Aster-only' subimage adjacent to Southern edge of 60 degree bounls -s dary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #   - 'Aster-only' subimage adjacent to Northern edge of 60 degree boundary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
          inMosaic      <- raster("./AsterSrtmBoth3ArcSecSub.tif")
     # Use these objects to get the extents of the subsets within the mosaic
          inAsterNorth  <- raster("./AsterBdyTestAbove60Sub.tif")
          inAsterSouth  <- raster("./AsterBdyTestBelow60Sub.tif")
     # Better, read them into a rasterStack
          mosaicLayers <- stack("./AsterSRTMBothTestSub.tif",
                                "./AsterBdyTestFGAbove60Sub.tif",
                                "./AsterBdyTestFGBelow60Sub.tif")
+    #
     # Create copy of input raster that we will use to create a 'fixed boundary' iamge
     # Even SIMPLER, according to raster documentation
+    #
          outMosaic <- raster(mosaicLayers[[1]])
          sFixedRasterFile <- "TestOutputRasterBdyFix.tif"
     # First, get the extent of the 'below 60' sub image in the (big) input raster
          northAsterEx <- extent(inAsterNorth)
          southAsterEx <- extent(inAsterSouth)
     # Get the values from the mosaic for this extent
          southCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],southAsterEx)
          northCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],northAsterEx)
     # Within the large input raster, we need the index of the first row of the 'Below 60 Aster' subimage.
     # Our plan: to replace this portion of the input mosaic with a linear combination of the original
     # input mosaic COPY (outMosaic) and the 'below 60' raster.
          firstNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[1])
          lastNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[(length(northCellsOfInterest) - 1)])
          northNrowsToProcess <- nrow(inAsterNorth)
          firstSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[1])
          lastSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[(length(southCellsOfInterest) - 1)])
          southNrowsToProcess <- nrow(inAsterSouth)
     #  create the output raster by copying the input file with 'block copy'
     message("Create output (fixed) mosaic")
     #browser()
          bs <- blockSize(mosaicLayers[[1]])
          outMosaic <- writeStart(outMosaic,sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
          for (iCtr in 1 : bs$n)
+         {
             mosaicVals <- getValues(mosaicLayers[[1]],row=bs$row[iCtr],nrows=bs$nrows[iCtr])
             writeValues(outMosaic,mosaicVals,bs$row[iCtr])
             message(sprintf(".....processed input mosaic block %d",iCtr))
+         }
          outMosaic <- writeStop(outMosaic)
     message("Input mosaic copied to output raster: now, process boundary")
     #browser()
     # now, we SHOULD be able to modify outMosaic with new 'column values'.
     # note: last north image and first south image rows are the same.
          northBorderEdgeRow <- lastNorthSubImgRow
          southBorderEdgeRow <- firstSouthSubImgRow + 1
     # The border 'edge' row is one row below the top of the south sub image
          southBrdrRowVals <- getValues(outMosaic,southBorderEdgeRow,1)
          northBrdrRowVals <- getValues(outMosaic,northBorderEdgeRow,1)
          brdrRowEdgeDelta <- southBrdrRowVals - northBrdrRowVals
     # Process the mosaic 'column-by-column'
     # For each column, extract the LAST 'numBlendRows cells adjacent to the bottomi
     # (southernmost) edge of the ASTER component, next to the SRTM border.
     # Then, moving 'north' along the colum of values, add the current colunn's 'ledge'
     # offset, attenuated by the exponential decay factor for that row.
     # Thus, the SRTM/ASTER 'ledge' is attenuated, the impact is greatest at the boundary,
     # least at the threshold, when the decay factor (empirically determined) goes to zero.
+    #
     # Add Exponential Decay constant: empirically selected (try a few, ee if you like result
     # Maybe plot it, see how many rows it takes for the 'effect' of SRTM goes to zero.
          dKc <- .045 # empirically determined constant, so that exp - 1 decay goes to zero in the desired number of image rows.
          numBlendRows <- 60
          begBlendRow <- (northBorderEdgeRow - numBlendRows) + 1
     #     srows <- seq(southBorderEdgeRow,((southBorderEdgeRow + numBlendRows) - 1))
          srows <- seq((begBlendRow),northBorderEdgeRow)
+    #
          for (curMosCol in 1 : ncol(mosaicLayers[[1]]))
+         {
     message(sprintf("transforming cur col: %d",curMosCol))
     #browser()
     #        deltaInc <- brdrRowEdgeDelta[curMosCol] / numBlendRows
             colVecCells <- cellFromRowCol(outMosaic,srows,curMosCol)
             colVecValues <- getValuesBlock(outMosaic,row=begBlendRow,nrows=numBlendRows,col=curMosCol,ncols=1)
             compareVec <- getValuesBlock(outMosaic,row=begBlendRow,nrows=numBlendRows+2,col=curMosCol,ncols=1)
             curRowBoundaryOffset <- brdrRowEdgeDelta[curMosCol]
     # note: we need to deal with NA values, which are possible.
     #        sumDeltaInc <- 0
     #        for (iCtr in numBlendRows : 1) # this moves from 'south to north', decay increases with distance.
             for (iCtr in 1 : numBlendRows) # this moves from 'south to north', decay increases with distance.
     #        for (iCtr in length(colVecValues) : 1) # this moves from 'south to north', decay increases with distance.
+            {
     # The idea: the closer to the border, the larger an increment we assign.
     # in any case, increment the offset so that it is correct for any iupcoming non-NA column value
                 if (!is.na(colVecValues[1]))
+                {
     #message(sprintf("cur col: %d Found a NON-NA vector",curMosCol))
     #browser()
     # exponential decay diminishes impact of the 'delta' edge adjustment with distance from the boundary.
                    dK <- exp(-(iCtr * dKc))
                    colToMod <- (numBlendRows-iCtr)+1
     #               qq <- as.integer(round(colVecValues[colToMod] + (curRowBoundaryOffset * dK)))
                    colVecValues[colToMod] <- as.integer(round(colVecValues[colToMod] + (curRowBoundaryOffset * dK)))
     #               colVecValues[iCtr] <- as.integer(round(colVecValues[iCtr] + (curRowBoundaryOffset * dK)))
+                }
+            }
     # Insert this vector back into the mosaic: this technique adopted from raster PDF doc 'replacement' help.
     # Looks like: I have to specify the row index (vector) to insert the values into the mosaic.
     message("write columns to out mosaic")
     #rowser()
             outMosaic[colVecCells] <- colVecValues
+         }
     # write the modified outMosaic values to the image file that we created.
     message("Done with raster blending : hit key to update raster disk file..")
     browser()
          writeRaster(outMosaic, sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
+    }

     #####################################################################################
     # DEMBoundarySolutionsBlock - Block version, uses a BLOCK write to create output mosaic
+    #
     # Collection of functions that attempt to mitigate ASTER / SRTM the boundary edge problem
+    #
+    #
+    #
     # Inputs:
     #   1) Blended Aster/SRTM DEM mosaic image.
     #   2) 'Aster-Only' subimage corresponding to the NORTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the SRTM component of the mosaic.
     #   3) 'Aster-Only' subimage corresponding to the SOUTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the ASTER component of the mosaic.
     #   4) Other parameters, as they are needed
+    #
     # To run:
+    #
     #  1) start R
     #  2) > source("DEMBoundarySolutions.r")
     #  3) > BlendAsterDEMPixelsBlock()
+    #
     # Author: Rick Reeves, NCEAS
     # May 31, 2011
     ##############################################################################################
     BlendAsterDEMPixelsBlock <- function()
+    {
       require(raster)
       require(rgdal)
     # Read the mosaic image and component sub-images:
     #   - Aster/SRTM mosaic image (with boundary edge to be fixed)
     #   - 'Aster-only' subimage adjacent to Southern edge of 60 degree bounls -s dary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #   - 'Aster-only' subimage adjacent to Northern edge of 60 degree boundary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #     inMosaic      <- raster("t:/rcr/BoundaryFixTest/AsterSrtmBoth3ArcSecSub.tif")
     #     inAsterNorth  <- raster("t:/rcr/BoundaryFixTest/AsterBdyTest3ArcSecFullGridAbove60.tif")
     #    inAsterSouth  <- raster("t:/rcr/BoundaryFixTest/AsterBdyTest3ArcSecFullGridBelow60.tif")
          inMosaic      <- raster("./AsterSrtmBoth3ArcSecSub.tif")
     # Use these objects to get the extents of the subsets within the mosaic
          inAsterNorth  <- raster("./AsterBdyTestAbove60Sub.tif")
          inAsterSouth  <- raster("./AsterBdyTestBelow60Sub.tif")
     # Better, read them into a rasterStack
     #     mosaicLayers <- stack("./AsterSrtmBoth3ArcSecSub.tif",
     #                           "./AsterBdyTest3ArcSecFullGridAbove60.tif",
     #                           "./AsterBdyTest3ArcSecFullGridBelow60.tif") # fix by exporting from ArcGIS
          mosaicLayers <- stack("./AsterSRTMBothTestSub.tif",
                                "./AsterBdyTestFGAbove60Sub.tif",
                                "./AsterBdyTestFGBelow60Sub.tif")
+    #
     # Create copy of input raster that we will use to create a 'fixed boundary' iamge
     # Even SIMPLER, according to raster documentation
+    #
          outMosaic <- raster(mosaicLayers[[1]])
          sFixedRasterFile <- "TestOutputRasterBdyFix.tif"
     # First, get the extent of the 'below 60' sub image in the (big) input raster
          northAsterEx <- extent(inAsterNorth)
          southAsterEx <- extent(inAsterSouth)
     # Get the values from the mosaic for this extent
          southCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],southAsterEx)
          northCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],northAsterEx)
     # Within the large input raster, we need the index of the first row of the 'Below 60 Aster' subimage.
     # Our plan: to replace this portion of the input mosaic with a linear combination of the original
     # input mosaic COPY (outMosaic) and the 'below 60' raster.
          firstNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[1])
          lastNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[(length(northCellsOfInterest) - 1)])
          northNrowsToProcess <- nrow(inAsterNorth)
          firstSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[1])
          lastSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[(length(southCellsOfInterest) - 1)])
          southNrowsToProcess <- nrow(inAsterSouth)
     #  create the output raster by copying the input file with 'block copy'
     message("Create output (fixed) mosaic")
     #browser()
          bs <- blockSize(mosaicLayers[[1]])
          outMosaic <- writeStart(outMosaic,sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
          for (iCtr in 1 : bs$n)
+         {
             mosaicVals <- getValues(mosaicLayers[[1]],row=bs$row[iCtr],nrows=bs$nrows[iCtr])
             writeValues(outMosaic,mosaicVals,bs$row[iCtr])
             message(sprintf(".....processed input mosaic block %d",iCtr))
+         }
          outMosaic <- writeStop(outMosaic)
     message("Input mosaic copied to output raster: now, process boundary")
     #browser()
     # now, we SHOULD be able to modify outMosaic with new 'column values'.
     # note: last north image and first south image rows are the same.
          northBorderEdgeRow <- lastNorthSubImgRow
          southBorderEdgeRow <- firstSouthSubImgRow + 1
     # The border 'edge' row is one row below the top of the south sub image
          southBrdrRowVals <- getValues(outMosaic,southBorderEdgeRow,1)
          northBrdrRowVals <- getValues(outMosaic,northBorderEdgeRow,1)
          brdrRowEdgeDelta <- southBrdrRowVals - northBrdrRowVals
     # Process the mosaic 'column-by-column'
     # For first effort, apply straightforward linear ramp function to the northernmost SRTM image rows.
     # In next iteration, we will modify the LAST 'numBlendRows' rows NORTH of the boundary with an 'extinction' function.
+    #
          numBlendRows <- 100
          srows <- seq(southBorderEdgeRow,((southBorderEdgeRow + numBlendRows) - 1))
+    #
          for (curMosCol in 1 : ncol(mosaicLayers[[1]]))
+         {
     message(sprintf("transforming cur col: %d",curMosCol))
     #browser()
             deltaInc <- brdrRowEdgeDelta[curMosCol] / numBlendRows
             colVecCells <- cellFromRowCol(outMosaic,srows,curMosCol)
             colVecValues <- getValuesBlock(outMosaic,row=southBorderEdgeRow,nrows=numBlendRows,col=curMosCol,ncols=1)
     # note: we need to deal with NA values, which are possible.
             sumDeltaInc <- 0
             for (iCtr in length(colVecValues) : 1)
+            {
     # The idea: the closer to the border, the larger an increment we assign.
     # in any case, increment the offset so that it is correct for any iupcoming non-NA column value
                 sumDeltaInc <- sumDeltaInc + deltaInc
                 if (!is.na(colVecValues[1]))
+                {
     #message(sprintf("cur col: %d Found a NON-NA vector",curMosCol))
     #browser()
                    colVecValues[iCtr] <- as.integer(round(colVecValues[iCtr] - sumDeltaInc))
+                }
+            }
     # Insert this vector back into the mosaic: this technique adopted from raster PDF doc 'replacement' help.
     # Looks like: I have to specify the row index (vector) to insert the values into the mosaic.
     #message("write columns to out mosaic")
     #browser()
             outMosaic[colVecCells] <- colVecValues
     #        outMosaic[srows,curMosCol] <- colVecValues #  this did not seem to work - june 1
+         }
     #message(sprintf("cur col: %d about to write to outMosaic",curMosCol))
     #browser()
     # write the modified outMosaic values to the image file that we created.
     message("Done with raster blending : hit key to update raster disk file..")
     browser()
          writeRaster(outMosaic, sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
+    }

     #####################################################################################
     # FixDEMMosaicBlendedRows: Implements Image Blending treatment to DEM Mosaic image
     # Collection of functions that attempt to mitigate ASTER / SRTM the boundary edge artifact
+    #
     # Inputs:
     #   1) DEM Mosaic image containing 'boundary edge' artifact
     #   2) 'Aster-Only' subimage corresponding to the NORTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the SRTM component of the mosaic.
     #   3) 'Aster-Only' subimage corresponding to the SOUTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the ASTER component of the mosaic.
     #   4) numBlendRows: number of number of rows south of 60 degrees N Lat to modify ('blend') by
     #      generating (lnear) combination of SRTM and ASTER image pixels.
+    #
     # Output:
     #   1) A copy of the input DEM mosaic, with 'blended image' boundary edge treatment.
+    #
     # To run:
+    #
     #  1) start R
     #  2) > source("FixDEMMosaicBlendedRows.r")
     #  3) > FixDEMMosaicBlendedRows()
+    #
     # Note: the input images used by this program located on /jupiter: /data/project/organisms/rcr/BoundaryFixTest
+    #
     # Author: Rick Reeves, NCEAS
     # June 12, 2011
+    #
     ##############################################################################################
     FixDEMMosaicBlendedRows <- function()
+    {
       require(raster)
       require(rgdal)
     # Key parameter: number of rows south of 60 degrees N Lat to modify ('blend') by
     # generating (lnear) combination of SRTM and ASTER image pixels.
       numBlendRows <- 100
     # Read the mosaic image and component sub-images:
     #   - Aster/SRTM mosaic image (with boundary edge to be fixed)
     #   - 'Aster-only' subimage adjacent to Southern edge of 60 degree bounls -s dary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #   - 'Aster-only' subimage adjacent to Northern edge of 60 degree boundary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #     inAsterNorthFG  <- raster("./aster60DegNorthRasterFG.tif")
     #     inAsterSouthFG  <- raster("./aster59DegSouthRasterFG.tif")
     #     inBdyAsterBoth  <- raster("./boundaryTestCase115DegBothAsterMos.img")
          inMosaic      <- raster("./AsterSrtm3ArcSecTestImg.tif")
          inAsterSouth  <- raster("./boundaryTest115_59DegreeAster.img")
          inAsterSouthFG  <- raster("./aster59DegSouthRasterFG.tif")
          inAsterNorth  <- raster("./boundaryTest115_60DegreeAster.img")
     # First, get the extent of the 'below 60' sub image in the (big) input raster
          northAsterEx <- extent(inAsterNorth)
          southAsterEx <- extent(inAsterSouth)
     # Create copy of input raster that we will use to create a 'fixed boundary' iamge
     # Even SIMPLER, according to raster documentation
          outMosaic <- raster(inMosaic)
          sFixedRasterFile <- "TestOutRasterBdyBlendFixAft2.tif"
     #  create the output raster by copying the input file with 'block copy'
     message("Copy input mosaic: Create output (repaired) mosaic")
     #browser()
          bs <- blockSize(inMosaic)
          outMosaic <- writeStart(outMosaic,sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
          for (iCtr in 1 : bs$n)
+         {
             mosaicVals <- getValues(inMosaic,row=bs$row[iCtr],nrows=bs$nrows[iCtr])
             writeValues(outMosaic,mosaicVals,bs$row[iCtr])
             message(sprintf(".....processed input mosaic block %d",iCtr))
+         }
          outMosaic <- writeStop(outMosaic)
+    #
     message("Input mosaic copied to output raster: now, process boundary")
     browser()
     # now, we SHOULD be able to modify outMosaic with new 'column values'.
          southCellsOfInterest <- cellsFromExtent(outMosaic,southAsterEx)
          northCellsOfInterest <- cellsFromExtent(outMosaic,northAsterEx)
          firstNorthSubImgRow <- rowFromCell(outMosaic,northCellsOfInterest[1])
          lastNorthSubImgRow <- rowFromCell(outMosaic,northCellsOfInterest[(length(northCellsOfInterest) - 1)])
          firstSouthSubImgRow <- rowFromCell(outMosaic,southCellsOfInterest[1])
          lastSouthSubImgRow <- rowFromCell(outMosaic,southCellsOfInterest[(length(southCellsOfInterest) - 1)])
     # note: last north image and first south image rows are the same.
          northBorderEdgeRow <- lastNorthSubImgRow
          southBorderEdgeRow <- firstSouthSubImgRow + 1
     # The border 'edge' row is one row below the top of the south sub image
          northBrdrRowVals <- getValues(outMosaic,northBorderEdgeRow,1)
          southBrdrRowVals <- getValues(outMosaic,southBorderEdgeRow,1)
     # Process the mosaic rows:
     # Blend the first 'numBlendRows' south of the 60 degree line (first rows of the SRTM component)
     # with the underlying and corresponding ASTER image rows.
     # For first attempt, use a linear combination; for second attempt, use an exponential function.
+    #
          nImageCols <- ncol(inMosaic)
          sRows <- seq(southBorderEdgeRow,((southBorderEdgeRow + numBlendRows) - 1))
+    #
          rowVecValuesBlend <- vector(mode="integer",length=nImageCols)
          iRowCtr <- 1
          for (curMosRow in sRows)
+         {
     message(sprintf("transforming cur row: %d",curMosRow))
     #:browser()
     # even simpler here, as we implement a 'linear distance decay' function,
     # and let each new column be a linear combination of the Aster and CGIAR image layers
             deltaInc <- iRowCtr / as.numeric(numBlendRows)
             colVecCells <- cellFromRowCol(outMosaic,curMosRow,1:nImageCols)
     # get the current row from the mosaic(SRTM data) and raster images
     # We get the ASTER values from a version of the 'South' ASTER image
     # with same extent as the image mosaic.
             rowVecValuesMosaic <- getValuesBlock(outMosaic,row=curMosRow,nrows=1,col=1,ncols=nImageCols)
             rowVecValuesAster <- getValuesBlock(inAsterSouthFG,row=curMosRow,nrows=1,col=1,ncols=nImageCols)
             rowVecValuesBlend <- (rowVecValuesMosaic * deltaInc) + (rowVecValuesAster * (1.0 - deltaInc))
     # create the merged row as a linear combination of the Mosaic and Aster image.
     # note: we need to deal with NA values, which are possible.
     # write the blended row to the output mosaic
             outMosaic[colVecCells] <- as.integer(round(rowVecValuesBlend))
             iRowCtr <- iRowCtr + 1
+         }
     #message(sprintf("cur col: %d about to write to outMosaic",curMosCol))
     #browser()
     # write the modified outMosaic values to the image file that we created.
     message("Done with raster blending : hit key to update raster disk file..")
     browser()
          writeRaster(outMosaic, sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
     #browser()
+    }

     #####################################################################################
     # FixDEMMosaicBlendedRowsExpDecay: Implements Image Blending treatment to DEM Mosaic image
     # Collection of functions that attempt to mitigate ASTER / SRTM the boundary edge artifact
     # This version implements the Exponential Decay function to control image blending rate.
     # Inputs:
     #   1) DEM Mosaic image containing 'boundary edge' artifact
     #   2) 'Aster-Only' subimage corresponding to the NORTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the SRTM component of the mosaic.
     #   3) 'Aster-Only' subimage corresponding to the SOUTHernmost part (adjacent to the 60 degree Lat
     #     boundary) of the ASTER component of the mosaic.
     #   4) numBlendRows: number of number of rows south of 60 degrees N Lat to modify ('blend') by
     #      generating (lnear) combination of SRTM and ASTER image pixels.
     #   5) dKc: Empirically-determined exponential decay constant. Rule of thumb: set this such
     #      that exponential decay function approaches zero as row 'numBlendRows' is treated.
+    #
     # Output:
     #   1) A copy of the input DEM mosaic, with 'blended image' boundary edge treatment.
+    #
     # To run:
+    #
     #  1) start R
     #  2) > source("FixDEMMosaicBlendedRows.r")
     #  3) > FixDEMMosaicBlendedRows()
+    #
     # Note: the input images used by this program located on /jupiter: /data/project/organisms/rcr/BoundaryFixTest
+    #
     # Author: Rick Reeves, NCEAS
     # June 13, 2011
+    #
     ##############################################################################################
     FixDEMMosaicBlendedRowsExpDecay <- function()
+    {
       require(raster)
       require(rgdal)
     # Key parameter: number of rows south of 60 degrees N Lat to modify ('blend') by
     # generating (lnear) combination of SRTM and ASTER image pixels.
       numBlendRows <- 100
       dKc <- .045 # empirically determined exponential decay constant. Select it
     #              so that exp - 1 decay goes to zero in the desired number of image rows.
     # Read the mosaic image and component sub-images:
     #   - Aster/SRTM mosaic image (with boundary edge to be fixed)
     #   - 'Aster-only' subimage adjacent to Southern edge of 60 degree bounls -s dary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #   - 'Aster-only' subimage adjacent to Northern edge of 60 degree boundary,
     #     sampled into a grid with same extent and spatial resolution as image mosaic
     #     inAsterNorthFG  <- raster("./aster60DegNorthRasterFG.tif")
     #     inAsterSouthFG  <- raster("./aster59DegSouthRasterFG.tif")
     #     inBdyAsterBoth  <- raster("./boundaryTestCase115DegBothAsterMos.img")
          inMosaic      <- raster("./AsterSrtm3ArcSecTestImg.tif")
          inAsterSouth  <- raster("./boundaryTest115_59DegreeAster.img")
          inAsterSouthFG  <- raster("./aster59DegSouthRasterFG.tif")
          inAsterNorth  <- raster("./boundaryTest115_60DegreeAster.img")
     # First, get the extent of the 'below 60' sub image in the (big) input raster
          northAsterEx <- extent(inAsterNorth)
          southAsterEx <- extent(inAsterSouth)
     # Create copy of input raster that we will use to create a 'fixed boundary' iamge
     # Even SIMPLER, according to raster documentation
          outMosaic <- raster(inMosaic)
          sFixedRasterFile <- "TestOutRasterBdyBlendFixExpDec.tif"
     #  create the output raster by copying the input file with 'block copy'
     message("Copy input mosaic: Create output (repaired) mosaic")
     #browser()
          bs <- blockSize(inMosaic)
          outMosaic <- writeStart(outMosaic,sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
          for (iCtr in 1 : bs$n)
+         {
             mosaicVals <- getValues(inMosaic,row=bs$row[iCtr],nrows=bs$nrows[iCtr])
             writeValues(outMosaic,mosaicVals,bs$row[iCtr])
             message(sprintf(".....processed input mosaic block %d",iCtr))
+         }
          outMosaic <- writeStop(outMosaic)
+    #
     message("Input mosaic copied to output raster: now, process boundary")
     browser()
     # now, we SHOULD be able to modify outMosaic with new 'column values'.
          southCellsOfInterest <- cellsFromExtent(outMosaic,southAsterEx)
          northCellsOfInterest <- cellsFromExtent(outMosaic,northAsterEx)
          firstNorthSubImgRow <- rowFromCell(outMosaic,northCellsOfInterest[1])
          lastNorthSubImgRow <- rowFromCell(outMosaic,northCellsOfInterest[(length(northCellsOfInterest) - 1)])
          firstSouthSubImgRow <- rowFromCell(outMosaic,southCellsOfInterest[1])
          lastSouthSubImgRow <- rowFromCell(outMosaic,southCellsOfInterest[(length(southCellsOfInterest) - 1)])
     # note: last north image and first south image rows are the same.
          northBorderEdgeRow <- lastNorthSubImgRow
          southBorderEdgeRow <- firstSouthSubImgRow + 1
     # The border 'edge' row is one row below the top of the south sub image
          northBrdrRowVals <- getValues(outMosaic,northBorderEdgeRow,1)
          southBrdrRowVals <- getValues(outMosaic,southBorderEdgeRow,1)
     # Process the mosaic rows:
     # Blend the first 'numBlendRows' south of the 60 degree line (first rows of the SRTM component)
     # with the underlying and corresponding ASTER image rows.
     # For first attempt, use a linear combination; for second attempt, use an exponential function.
+    #
          nImageCols <- ncol(inMosaic)
          sRows <- seq(southBorderEdgeRow,((southBorderEdgeRow + numBlendRows) - 1))
+    #
          rowVecValuesBlend <- vector(mode="integer",length=nImageCols)
          iRowCtr <- 1
          for (curMosRow in sRows)
+         {
              message(sprintf("transforming cur row: %d",curMosRow))
     # exponential decay diminishes impact of the 'delta' edge adjustment with distance from the boundary.
              dK <- exp(-(iRowCtr * dKc))
     #         colToMod <- (numBlendRows-iCtr)+1
     #               qq <- as.integer(round(colVecValues[colToMod] + (curRowBoundaryOffset * dK)))
     #               colVecValues[colToMod] <- as.integer(round(colVecValues[colToMod] + (curRowBoundaryOffset * dK)))
     #               colVecValues[iCtr] <- as.integer(round(colVecValues[iCtr] + (curRowBoundaryOffset * dK)))

... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff

Project

General

Profile

Revision a09a09c0

Added by Jim Regetz over 12 years ago