Revision 80390607
Added by Rick Reeves over 13 years ago
- ID 80390607717cb5b1e4a951a017b49c7163cb4bc4
| terrain/rscripts/MosaicCdemDifferencePlots.r | ||
|---|---|---|
| 1 |
#setwd("f:/Projects/NCEAS/EnviromentAndOrganisms")
|
|
| 2 |
CreateDiffPlots <- function() |
|
| 3 |
{
|
|
| 1 |
############################################################################################## |
|
| 2 |
# |
|
| 3 |
# CreateImageDiffPlots |
|
| 4 |
# |
|
| 5 |
# R script generates collection plots that explore elevation discontinuities |
|
| 6 |
# border between ASTER and STRM data in mosaic images under development for Iplant. |
|
| 7 |
# Plots overwritten to single graphics device for first version. |
|
| 8 |
# Panel Plots to be added to next version. |
|
| 9 |
# |
|
| 10 |
# To run: |
|
| 11 |
# |
|
| 12 |
|
|
| 13 |
# 1) start R, |
|
| 14 |
# |
|
| 15 |
# |
|
| 16 |
# Input file: Csv table containing elevations for randomly-selected point pairs from |
|
| 17 |
# ASTER/CGIAR mosaic and CDEM images. Table is generated by the script: 'makeImagePairTable.r' |
|
| 18 |
# Author: Rick Reeves, NCEAS |
|
| 19 |
# May 1, 2011 |
|
| 20 |
############################################################################################## |
|
| 21 |
#CreateImageDiffPlots <- function() |
|
| 22 |
#{
|
|
| 23 |
|
|
| 24 |
# Read input table that was sorted OFFLINE in Excel on the first column (ColumnID) |
|
| 25 |
|
|
| 4 | 26 |
pointTable <-read.csv("tableForMark4000_5_8_SortColID.csv")
|
| 5 | 27 |
# |
| 6 | 28 |
# Table created by randomly sampling from two superimposed |
| ... | ... | |
| 49 | 71 |
# Difference between Mosaic (ASTER or CGIAR or border) |
| 50 | 72 |
# and CDEM elevation as pertentage of the mosaic elevation |
| 51 | 73 |
|
| 52 |
pointTable <-cbind(pointTable,pointTable$diffNorthMosaicCDEM/pointTable$elevNorth) * 100
|
|
| 53 |
pointTable <-cbind(pointTable,pointTable$diffSouthMosaicCDEM/pointTable$elevSouth) * 100
|
|
| 74 |
pointTable <-cbind(pointTable,(pointTable$diffNorthMosaicCDEM/pointTable$elevNorth * 100))
|
|
| 75 |
pointTable <-cbind(pointTable,(pointTable$diffSouthMosaicCDEM/pointTable$elevSouth * 100))
|
|
| 54 | 76 |
|
| 55 | 77 |
colnames(pointTable)[10] <- "magDiffMosaicCDEMNorthPct" |
| 56 | 78 |
colnames(pointTable)[11] <- "magDiffMosaicCDEMSouthPct" |
| ... | ... | |
| 63 | 85 |
northRowTbl = pointTable[northRowIndexes,] |
| 64 | 86 |
borderRowTbl = pointTable[borderRowIndexes,] |
| 65 | 87 |
southRowTbl = pointTable[southRowIndexes,] |
| 66 |
message("hit key to create plots...")
|
|
| 67 |
browser() |
|
| 88 |
|
|
| 68 | 89 |
# |
| 69 | 90 |
# Create plots |
| 70 | 91 |
# |
| 71 | 92 |
# First, plot North/South pixel elevation for North/Border/South subsets 'distance East from Western edge'. |
| 72 | 93 |
# |
| 73 | 94 |
plot(northRowTbl$ColumnID,northRowTbl$elevNorth,main="North (ASTER) Mosaic North Pixel Elevation",col="red",pch=".") |
| 95 |
message("Hit key for next plot...")
|
|
| 96 |
browser() |
|
| 74 | 97 |
plot(borderRowTbl$ColumnID,borderRowTbl$elevNorth,main="Border Mosaic North Elevation",col="darkgreen",pch=".") |
| 98 |
message("Hit key for next plot...")
|
|
| 99 |
browser() |
|
| 75 | 100 |
plot(southRowTbl$ColumnID,southRowTbl$elevNorth,main="South (SRTM) Mosaic North Elevation",col="blue",pch=".") |
| 76 |
|
|
| 101 |
message("Hit key for next plot...")
|
|
| 102 |
browser() |
|
| 77 | 103 |
plot(northRowTbl$ColumnID,northRowTbl$elevNorth,main="North (ASTER) CDEM North Pixel Elevation",col="red",pch=".") |
| 104 |
message("Hit key for next plot...")
|
|
| 105 |
browser() |
|
| 78 | 106 |
plot(borderRowTbl$ColumnID,borderRowTbl$elevNorth,main="Border CDEM North Pixel Elevation",col="darkgreen",pch=".") |
| 107 |
message("Hit key for next plot...")
|
|
| 108 |
browser() |
|
| 79 | 109 |
plot(southRowTbl$ColumnID,southRowTbl$elevNorth,main="South (SRTM) CDEM North Elevation",col="blue",pch=".") |
| 80 |
|
|
| 110 |
message("Hit key for next plot...")
|
|
| 111 |
browser() |
|
| 81 | 112 |
plot(northRowTbl$ColumnID,northRowTbl$elevSouth,main="North (ASTER) Mosaic South Pixel Elevation",col="red",pch=".") |
| 113 |
message("Hit key for next plot...")
|
|
| 114 |
browser() |
|
| 82 | 115 |
plot(borderRowTbl$ColumnID,borderRowTbl$elevSouth,main="Border Mosaic South Elevation",col="darkgreen",pch=".") |
| 116 |
message("Hit key for next plot...")
|
|
| 117 |
browser() |
|
| 83 | 118 |
plot(southRowTbl$ColumnID,southRowTbl$elevSouth,main="South (SRTM) Mosaic South Elevation",col="blue",pch=".") |
| 84 |
|
|
| 119 |
message("Hit key for next plot...")
|
|
| 120 |
browser() |
|
| 85 | 121 |
plot(northRowTbl$ColumnID,northRowTbl$elevSouth,main="North (ASTER) CDEM South Pixel Elevation",col="red",pch=".") |
| 122 |
message("Hit key for next plot...")
|
|
| 123 |
browser() |
|
| 86 | 124 |
plot(borderRowTbl$ColumnID,borderRowTbl$elevSouth,main="Border CDEM South Pixel Elevation",col="darkgreen",pch=".") |
| 125 |
message("Hit key for next plot...")
|
|
| 126 |
browser() |
|
| 87 | 127 |
plot(southRowTbl$ColumnID,southRowTbl$elevSouth,main="South (SRTM) CDEM South Elevation",col="blue",pch=".") |
| 128 |
message("Hit key for next plot...")
|
|
| 129 |
browser() |
|
| 88 | 130 |
# |
| 89 | 131 |
# Next, (for now) difference between North and South pixels in pair, for various regions |
| 90 | 132 |
# |
| 91 | 133 |
plot(northRowTbl$ColumnID,northRowTbl$diffMosaicNorthSouth,main="North (ASTER) Mosaic North/South Diff",col="red",pch=".") |
| 134 |
message("Hit key for next plot...")
|
|
| 135 |
browser() |
|
| 92 | 136 |
plot(borderRowTbl$ColumnID,borderRowTbl$diffMosaicNorthSouth,main="Border Mosaic North/South Diff",col="darkgreen",pch=".") |
| 137 |
message("Hit key for next plot...")
|
|
| 138 |
browser() |
|
| 93 | 139 |
plot(southRowTbl$ColumnID,southRowTbl$diffMosaicNorthSouth,main="South (SRTM) Mosaic North/South Diff",col="darkgreen",pch=".") |
| 140 |
message("Hit key for next plot...")
|
|
| 141 |
browser() |
|
| 94 | 142 |
# |
| 95 | 143 |
plot(northRowTbl$ColumnID,northRowTbl$diffCDEMNorthSouth,main="North (ASTER) CDEM North/South Diff",col="red",pch=".") |
| 144 |
message("Hit key for next plot...")
|
|
| 145 |
browser() |
|
| 96 | 146 |
plot(borderRowTbl$ColumnID,borderRowTbl$diffCDEMNorthSouth,main="Border CDEM North/South Diff",col="darkgreen",pch=".") |
| 147 |
message("Hit key for next plot...")
|
|
| 148 |
browser() |
|
| 97 | 149 |
plot(southRowTbl$ColumnID,southRowTbl$diffCDEMNorthSouth,main="South (SRTM) CDEM North/South Diff",col="blue",pch=".") |
| 150 |
message("Hit key for next plot...")
|
|
| 151 |
browser() |
|
| 98 | 152 |
|
| 99 | 153 |
# Next: the Mosaic/CDEM difference divided by the (ASTER or SRTM) elevation |
| 100 |
|
|
| 154 |
|
|
| 101 | 155 |
plot(northRowTbl$ColumnID,northRowTbl$magDiffMosaicCDEMNorthPct,col="red",pch=".",main="Elev Diff as PerCent of Mosaic Elev (North)") |
| 156 |
message("Hit key for next plot...")
|
|
| 157 |
browser() |
|
| 102 | 158 |
plot(borderRowTbl$ColumnID,borderRowTbl$magDiffMosaicCDEMNorthPct,col="darkgreen",pch=".",main="Elev Diff as PerCent of Mosaic Elev (Border)") |
| 159 |
message("Hit key for next plot...")
|
|
| 160 |
browser() |
|
| 103 | 161 |
plot(southRowTbl$ColumnID,southRowTbl$magDiffMosaicCDEMSouthPct,col="blue",pch=".",main="Elev Diff as PerCent of Mosaic Elev (South)") |
| 104 |
|
|
| 162 |
message("Hit key for next plot...")
|
|
| 163 |
browser() |
|
| 105 | 164 |
# Finally: Mosaic vs CDEM difference vs elevation (north or south) |
| 106 | 165 |
|
| 107 | 166 |
plot(northRowTbl$elevNorth,northRowTbl$diffCDEMNorthSouth,col="red",pch=".",main="ASTER Elev Diff vs Mosaic Elev (North Pixel)") |
| 108 | 167 |
# plot(northRowTbl$elevSouth,southRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="ASTER Elev Diff vs Mosaic Elev (South Pixel)") |
| 109 |
|
|
| 168 |
message("Hit key for next plot...")
|
|
| 169 |
browser() |
|
| 110 | 170 |
plot(borderRowTbl$elevNorth,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="Border Elev Diff vs Mosaic Elev (North Pixel)") |
| 111 | 171 |
# plot(borderRowTbl$elevSouth,borderRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="Border Elev Diff vs Mosaic Elev (South Pixel)") |
| 112 |
|
|
| 172 |
message("Hit key for next plot...")
|
|
| 173 |
browser() |
|
| 113 | 174 |
plot(southRowTbl$elevNorth,southRowTbl$diffCDEMNorthSouth,col="blue",pch=".",main="SRTM Elev Diff vs Mosaic Elev (North Pixel)") |
| 114 | 175 |
# plot(southRowTbl$elevSouth,southRowTbl$diffCDEMNorthSouth,col="darkgreen",pch=".",main="SRTM Elev Diff vs Mosaic Elev (South Pixel)") |
| 115 |
|
|
| 116 |
|
|
| 117 |
|
|
| 118 |
|
|
| 119 |
|
|
| 120 |
|
|
| 121 |
} |
|
| 176 |
message("end of plots")
|
|
| 177 |
#} |
|
Also available in: Unified diff
First working version - very simple - just 'source' the file and press <CR> to walk through the plots.
Next version: 'panel plots' and improved plot axes.