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

#

# extractDemSubimages

# 

# R script extracts and saves subimages from various inputs to the 'image boundary analysis';

# also produces 'difference images' by subtracting CDEM components from other DEM components.

# Output files written as GeoTiff files. 

#

# Author: Rick Reeves, NCEAS

# May 9, 2011

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

#

extractDemSubimages <- function()

{

require(raster)

require(rgdal)

#inputRasterMerge <- raster(sFirstImageName)

#inputSecondRaster <- raster(sSecondImageName)

inputRasterAster <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdyASTER_BLEven.tif")

inputRasterSRTM <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdySRTM_BLEven.tif")

inputRasterMerge <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")

inputRasterCDEM <- raster("/data/project/organisms/rcr/ValidateBoundary/CDEMMosCgiarAsterBdy_BLEven.tif")

inputCdemAster24Row <- raster("/data/project/organisms/rcr/ValidateBoundary/EdgeAsterCDEM24row.tif")

inputCdemSrtm24Row <- raster("/data/project/organisms/rcr/ValidateBoundary/EdgeSrtmCDEM24row.tif")

#

# Difference image for entire merged image takes a while to create, 

# so we created it once, now read it back in.

#

rDeltaWhole <- raster("/data/project/organisms/rcr/ValidateBoundary/DeltaEntireImage.tif")

rDeltaWhole@data@values <-getValues(rDeltaWhole)

# Create extent objects used to extract raster subimges. 

# The object will be centered along the 60 degree North latitude line,

# and have varying depths (number of rows).

# The Western Canada study area runs from -135 (west) to -100 (west) longitude,

# and 55.0 to 64.00 degrees (north) latitude. 

# the ASTER and SRTM/CGIAR image components are merged at the 60 Deg N Latitude line.

eTestAreaExtentAster <- extent(-135.0,-105.0, 60.00009,60.010) # Creates 12 row subimage North of border (All ASTER)

eTestAreaExtentBorder <- extent(-135.0,-105.0, 59.995,60.005) # Creates a 12 row subimage centered on border

eTestAreaExtentSRTM <- extent(-135.0,-105.0, 59.990,60.00) # Creates a 12 row subimage South of border (All SRTM)

eTestAreaExtentAster <- extent(-135.0,-105.0, 60.00009,60.020) # Creates 24 row subimage North of border (All ASTER)

eTestAreaExtentBorder <- extent(-135.0,-105.0, 59.99,60.010) # Creates a 24 row subimage centered on border

eTestAreaExtentSRTM <- extent(-135.0,-105.0, 59.98,60.00) # Creates a 24 row subimage South of border (All SRTM)

# Extract a sub image corresponding to the selected extent.

# Two different alternatives:

# The extract() function returns a vector of cell values, 

# the crop() function returns a complete raster* object.

message("extracts")

browser()

vEdgeRegionAster <- extract(inputRasterAster,eTestAreaExtentAster)

rEdgeRegionAster <- crop(inputRasterAster,eTestAreaExtentAster)

vEdgeRegionAsterDelta <- extract(rDeltaWhole,eTestAreaExtentAster)

rEdgeRegionAsterDelta <- crop(rDeltaWhole,eTestAreaExtentAster)

#vEdgeRegionAsterCDEM <- extract(inputRasterCDEM,eTestAreaExtentAster)

#rEdgeRegionAsterCDEM <- crop(inputRasterCDEM,eTestAreaExtentAster)

vEdgeRegionBorder <- extract(inputRasterMerge,eTestAreaExtentBorder)

rEdgeRegionBorder <- crop(inputRasterMerge,eTestAreaExtentBorder)

vEdgeRegionBorderDelta <- extract(rDeltaWhole,eTestAreaExtentBorder)

rEdgeRegionBorderDelta <- crop(rDeltaWhole,eTestAreaExtentBorder)

vEdgeRegionSRTM <- extract(inputRasterSRTM,eTestAreaExtentSRTM)

rEdgeRegionSRTM <- crop(inputRasterSRTM,eTestAreaExtentSRTM)

vEdgeRegionSRTMDelta <- extract(rDeltaWhole,eTestAreaExtentSRTM)

rEdgeRegionSRTMDelta <- crop(rDeltaWhole,eTestAreaExtentSRTM)

#vEdgeRegionSrtmCDEM <- extract(inputRasterCDEM,eTestAreaExtentSRTM)

#rEdgeRegionSrtmCDEM <- crop(inputRasterCDEM,eTestAreaExtentSRTM)

message("differences")

browser()

# 

# write these files to disk, get stats.

#

rDeltaEdgeAster <- rEdgeRegionAster - inputCdemAster24Row

rDeltaEdgeBorder <- rEdgeRegionBorder - rEdgeRegionBorderDelta

rDeltaEdgeSRTM <- rEdgeRegionSRTM - rEdgeRegionSrtmCDEM

# Important: In order for the image subtraction to work, the extents

#            of the two images must be IDENTICAL. I used ArcMap GIS Raster Crop By Mask

#            to create subimages with identical extents. 

# Computeexter the difference image  for the entire study area, and for the region along

# the boundary (narrow, maybe 10 pixels either side)

# rDeltaEdge <- rEdgeRegionFirst - rEdgeRegionSecond

# May 10 2011 note: One of these statements generates the 'different origin' fatal error. 

# my workaround was to create ALL of the subimages used here with gdalwarp, and then 

# read them in and subtract them. 

# see other R script: OnlyImage Differences.r

# I will create a reproducable example, send to Raster authors. 

rDeltaEdgeAster <- rEdgeRegionAster - inputCdemAster24Row

rDeltaEdgeSRTM <- rEdgeRegionSRTM - inputCdemSrtm24Row

#rDeltaWhole <- inputRasterMerge - inputRasterCDEM

#writeRaster(rDeltaWhole,filename="DeltaEntireImage.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)

# Using the large difference image, compute subimagee statistics for areas

# North (ASTER) and South (SRTM) of the boundary. These give us an idea 

# re: differences between ASTER and CDEM and CGIAR/SRT and CDEM

# what is raster package way of using subscripts to extract? 

# Now, the interesting part: using the boundary difference image, randomly select 

# one-degree N-S strips throughout the image, and compare adjacent pixel pairs 

# above and below the boundary with pixel pairs straddling the boundary. Subtract  

# the pairs, save the collection of (absolute value) of the differences in a vector,

# so that we have a population of differences above, below, and straddling the boundary 

# line. Compare the populations.

# get a vector of random column index numbers, constrained by column dimension of image

# Loop three times, sampling pixel pairs from above, below, across the border

nColsToGet <-20000

iDiffVecNorth <- vector(mode="integer",length=nColsToGet)

iDiffVecBorder <- vector(mode="integer",length=nColsToGet)

iDiffVecSouth <- vector(mode="integer",length=nColsToGet)

#colsToGet <-sample(1:50,nColsToGet)

# Note: initially, sample the same columns in all regions to get a profile.

#       other 'sample()' calls can be commented out to sample differenct

#       coluns in each 'region'.

# iDiffVecxxxx is a population of differences between adjacent cell pairs. 

# Compute iDiffVecNorth/Border/South on either side of border, and across it. 

# note that North and South samples taken from larger difference image for 

# entire mosaic (sub) image; iDiffBorder taken from the edge region extracted

# from the center of the lerger image.

# Remember, we are sampling a PAIR of pixels (same column from two adjacent rows)

# NOTR 5/9: need to rethink these sampling loops for this new version. Comment for now!

#colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)

#message("North Sample")

#browser()

# debug

#nColsToGet <- 2

#colsToGet <- c(20,100)

#iFirstRow <- 300

#iCtr = 1

#for (iNextCol in colsToGet)

#{

#  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)

#  neighborCells <- rDeltaWhole@data@values[rColVec]

#  iDiffVecNorth[iCtr] <- neighborCells[2] - neighborCells[1]

#  iCtr = iCtr + 1

#}

#

#message("Border Sample - different columns")

#browser()

#colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)

#iFirstRow <- 6 # straddle the border of 12 row center section 

#iCtr = 1

#for (iNextCol in colsToGet)

#{

#  rColVec <- cellFromRowCol(rDeltaEdge,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)

#  neighborCells <- rDeltaEdge@data@values[rColVec]

#  iDiffVecBorder[iCtr] <- neighborCells[2] - neighborCells[1]

#  iCtr = iCtr + 1 

#}

##

#message("South Sample - different columns")

##browser()

#colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)

#iFirstRow <- 3600

#iCtr = 1

#for (iNextCol in colsToGet)

#{

#  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)

#  neighborCells <- rDeltaWhole@data@values[rColVec]

#  iDiffVecSouth[iCtr] <- neighborCells[2] - neighborCells[1]

#  iCtr = iCtr + 1 

#}

## Compute iDiffVecs on either side of border, and across it. 

#message("Check the cell difference vectors...")

##browser()

# summary stats for each population

#sNorthSum <- sprintf("ASTER sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",

#                     min(iDiffVecNorth,na.rm=TRUE),median(iDiffVecNorth,na.rm=TRUE),mean(iDiffVecNorth,na.rm=TRUE),

#                     max(iDiffVecNorth,na.rm=TRUE),var(iDiffVecNorth,na.rm=TRUE),sd(iDiffVecNorth,na.rm=TRUE))

#

#sBorderSum <- sprintf("Border sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",

#                    min(iDiffVecBorder,na.rm=TRUE),median(iDiffVecBorder,na.rm=TRUE),mean(iDiffVecBorder,na.rm=TRUE),

#                    max(iDiffVecBorder,na.rm=TRUE),var(iDiffVecBorder,na.rm=TRUE),sd(iDiffVecBorder,na.rm=TRUE))

#sSouthSum <- sprintf("STRM sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",

#                     min(iDiffVecSouth,na.rm=TRUE),median(iDiffVecSouth,na.rm=TRUE),mean(iDiffVecSouth,na.rm=TRUE),

#                     max(iDiffVecSouth,na.rm=TRUE),var(iDiffVecSouth,na.rm=TRUE),sd(iDiffVecSouth,na.rm=TRUE))

##

#message(sprintf("statistics for %d N/S adjacent pixel pairs from three mosaic image regions:",nColsToGet))

#message(sNorthSum)

#message(sBorderSum)

#message(sSouthSum)

message("hit key to write output images...")

browser()

writeRaster(rDeltaEdgeSRTM,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionSRTMOnlyDelta24row.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)

writeRaster(rDeltaEdgeAster,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionAsterOnlyDelta24row.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)

# Create this image one time, read it in thereafter.

}

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

#

# R Script performs edge filtering on Raster files

#

#

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

filterDiffRasterImage <- function(sInCDEMRasterName,sInAstCgiarMosName)

{

require(raster)

require(rgdal)

# Read a raster image file into Raster object

#setwd("/data/project/organisms/rcr/ValidateBoundary")

#setwd("t:/organisms/rcr/ValidateBoundary")

#sInCDEMName <- "CanadaDemMosTuesdayNewWGS84.tif"

#sInAstCgiarMosName <- "mergeCgiarAsterBdyFinalBLMonday.tif"

inputCDEMRaster <- raster(sInCDEMRasterName)

inputAsterCgiarRaster <- raster(sInAstCgiarMosName)

inputAsterCgiarRasterRes <- raster("./mergeCgiarAsterTuesdayCDEMGrid2.tif")

# Create an extent object with the coordinates of the 'edge zone' 

# along the 60 degree North latitude line

# The Western Canada study area runs from -135 (west) to -100 (west) longitude,

# and 55.0 to 64.00 degrees (north) latitude. 

eTestAreaExtent <- extent(-135.2,-100.2, 59.995,60.005) # this results in a 12 row subimage

eTestAreaExtent <- extent(-135.2,-100.2, 59.997,60.001) # Creates a 5 row subimage CENTERED on 60N Latitude

# extract a sub area corresponding to the extent.

# The extract() function returns a vector of cell values, 

# the crop() function returns a raster* object.

vEdgeRegionCDEM <- extract(inputCDEMRaster,eTestAreaExtent)

rEdgeRegionCDEM <- crop(inputCDEMRaster,eTestAreaExtent)

vEdgeRegionAsterCgiar <- extract(inputAsterCgiarRaster,eTestAreaExtent)

rEdgeRegionAsterCgiar <- crop(inputAsterCgiarRaster,eTestAreaExtent)

vEdgeRegionAsterCgiarRes <- extract(inputAsterCgiarRasterRes,eTestAreaExtent)

rEdgeRegionAsterCgiarRes <- crop(inputAsterCgiarRasterRes,eTestAreaExtent)

# for raster subtraction to work, the raster object extracts must work. 

# for some reason, the extent for these two were off by .1 degree. 

# here is a workaround - It resulted in an 'aligned' image pair, for 

# which the difference image values matched the actual differences 

# between ASTER and CDEM values.

extent(rEdgeRegionAsterCgiar) <- extent(rEdgeRegionCDEM) # only do if extents dont match

rDelta <- rEdgeRegionAsterCgiar - rEdgeRegionCDEM

rDeltaRes <- rEdgeRegionAsterCgiarRes - rEdgeRegionCDEM

rBigDeltaRes <- inputAsterCgiarRasterRes - inputCDEMRaster

# compare the statistics between these two difference images

# Construct the filter kernal: a 3 x 3 integer matrix

# Execute the filter

# Compute image stats?

# display the file?

# Write the output file to disk

writeRaster(rEdgeRegionCDEM,filename="EdgeRegionCDEM5RowsS.tif",

            format="GTiff",datatype="INT2S",overwrite=TRUE)

writeRaster(rEdgeRegionAsterCgiar,filename="EdgeRegionAsterCgiar5RowsS.tif",

            format="GTiff",datatype="INT2S",overwrite=TRUE)

writeRaster(rEdgeRegionAsterCgiarRes,filename="EdgeRegionAsterCgiar5RowsResamp.tif",

format="GTiff",datatype="INT2S",overwrite=TRUE)

writeRaster(rDelta,filename="EdgeRegionDelta5RowsS.tif",

            format="GTiff",datatype="INT2S",overwrite=TRUE)

writeRaster(rDeltaRes,filename="EdgeRegionDelta5RowsResS.tif",

            format="GTiff",datatype="INT2S",overwrite=TRUE)

#writeRaster(rBigDelta,filename="CgiarCDEMRegionDelta.tif",

#            format="GTiff",datatype="INT2S",overwrite=TRUE)

#writeRaster(rBigDeltaRes,filename="CgiarCDEMRegionDeltaRes.tif",

#            format="GTiff",datatype="INT2S",overwrite=TRUE)            

}

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

#

# makeMarkTable.r

# 

# R script generates table of adjacent ASTER / SRTM / Mosaic Border / CDEM pixel values

# to be used to generate plots of 'delta elevation' vs 'pixel pair proximity' and 'elevation'

# suggested by Mark Schildhauer on May 3. 

#

# Author: Rick Reeves, NCEAS

# May 4, 2011

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

#

makeMarkTable <- function()

{

require(raster)

require(rgdal)

inputCgiarRaster  <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdySRTM_BL.tif")

inputAsterRaster  <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdyASTER_BL.tif")

inputMosaicRaster <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")

inputCDEMRaster   <- raster("/data/project/organisms/rcr/ValidateBoundary/CDemMosTuesdayClipMergeSpace.tif")

#

# Difference image for entire merged image takes a while to create, 

# so we created it once, now read it back in.

#

rDeltaWhole <- raster("/data/project/organisms/rcr/ValidateBoundary/DeltaEntireImage.tif")

rDeltaWhole@data@values <-getValues(rDeltaWhole)

# Create extent objects used to extract raster subimges. 

# The object will be centered along the 60 degree North latitude line,

# and have varying depths (number of rows).

# The Western Canada study area runs from -135 (west) to -100 (west) longitude,

# and 55.0 to 64.00 degrees (north) latitude. 

# the ASTER and SRTM/CGIAR image components are merged at the 60 Deg N Latitude line.

#eTestAreaExtent <- extent(-135.2,-100.2, 59.997,60.001) # Creates a 5 row subimage

eTestAreaExtent <- extent(-135.2,-100.2, 59.995,60.005) # Creates a 12 row subimage

# Extract a sub image corresponding to the selected extent.

# Two different alternatives:

# The extract() function returns a vector of cell values, 

# the crop() function returns a complete raster* object.

rEdgeRegionAster <- crop(inputAsterRaster,eTestAreaExtent)

rEdgeRegionCgiar <- crop(inputCgiarRaster,eTestAreaExtent)

rEdgeRegionMosaic <- crop(inputMosaicRaster,eTestAreaExtent)

rEdgeRegionCDEM <- crop(inputCDEMRaster,eTestAreaExtent)

# Important: In order for the image subtraction to work, the extents

#            of the two images must be IDENTICAL. I used ArcMap GIS Raster Crop By Mask

#            to create subimages with identical extents. 

# Compute the difference image  for the entire study area, and for the region along

# the boundary (narrow, maybe 10 pixels either side)

rEdgeRegionDelta <- rEdgeRegionMosaic - rEdgeRegionCDEM  # not used in this version (yet)

# get a vector of random column index numbers, constrained by column dimension of image

# Loop three times, sampling pixel pairs from above, below, across the border

nColsToGet <- 4000

iDiffVecNorth <- vector(mode="integer",length=nColsToGet)

iDiffVecBorder <- vector(mode="integer",length=nColsToGet)

iDiffVecSouth <- vector(mode="integer",length=nColsToGet)

# Here is the output matrix, (nColSamples * 3) rows, four columns 

colsToGet <-sample(1:inputMosaicRaster@ncols,nColsToGet)

mOutTable <- matrix(nrow=(nColsToGet * 3), ncol=5)

colnames(mOutTable) <- c("ColumnID","elevNorth","elevSouth","cdemNorth","cdemSouth")

# in this difference image, the border edge occurs at column 6

# NOTE Wed eve: A better start and end for the border: firstRow=5, lastRow = 8....

iFirstRow <- 5 #4 # two rows before the border           

iLastRow  <- 8 #7 # two rows after the border

iRowCtr = 1 # points to latest output table being written

message(sprintf("starting col sample loop"))

for (iNextCol in colsToGet)

{

#message(sprintf("in col sample loop getting col sample %d - hit key..",iNextCol))

#browser()

# slight misalignement in the Aster and Cgiar edge regions, getting NAs wheen extracting values. 

# For now, get values on both sides of border from the mosaic

#   rColVecAster  <- cellFromRowCol(rEdgeRegionAster,iFirstRow:(iLastRow),iNextCol:iNextCol) # figure out the actual row numbers for mosaic components LATER

#   rColVecCgiar  <- cellFromRowCol(rEdgeRegionCgiar,iFirstRow:(iLastRow),iNextCol:iNextCol)

#   rColVecAster  <- cellFromRowCol(rEdgeRegionMosaic,iFirstRow:(iLastRow),iNextCol:iNextCol) # approximation is the Cgiar/Aster area in mosaic.

#   rColVecCgiar  <- cellFromRowCol(rEdgeRegionMosaic,iFirstRow:(iLastRow),iNextCol:iNextCol)

   rColVecMosaic <- cellFromRowCol(rEdgeRegionMosaic,iFirstRow:(iLastRow),iNextCol:iNextCol)

   rColVecCDEM   <- cellFromRowCol(rEdgeRegionCDEM,iFirstRow:(iLastRow),iNextCol:iNextCol)  

# Split the column vector into the pairs that Mark requested 

# For each column sampled, the output table has three rows:

#

# Row 1: from (input) Aster layer:  Two pixels above border

# Row 2: from (input) Mosaic layer: Two pixels straddling border rDeltaWhole@data@values[rColVec]

# Row 3: from (input) Cgiar layer:  Two pixels below border

#

# each with four columns, arranged into two column pairs:

#

#  North Pixel and South Pixel elevation, North Pixel and South Pixel CDEM (baseline)

# 

   mOutTable[iRowCtr,1]   <- iNextCol           # The (radomly sampled) col ID (surrogate for longitude) 

   mOutTable[iRowCtr,2:3] <- rEdgeRegionMosaic@data@values[rColVecMosaic[1:2]]  # First column pair from extracted vector: 

   mOutTable[iRowCtr,4:5] <- rEdgeRegionCDEM@data@values[rColVecCDEM[1:2]]   # entirely top (ASTER part of mosaic) image 

   iRowCtr = iRowCtr + 1

#  

   mOutTable[iRowCtr,1]   <- iNextCol           # The (radomly sampled) col ID (surrogate for longitude)

   mOutTable[iRowCtr,2:3] <- rEdgeRegionMosaic@data@values[rColVecMosaic[2:3]] # Second column pair: 

   mOutTable[iRowCtr,4:5] <- rEdgeRegionCDEM@data@values[rColVecCDEM[2:3]]   # straddles border region 

   iRowCtr = iRowCtr + 1  

#   

   mOutTable[iRowCtr,1]   <- iNextCol           # The (radomly sampled) col ID (surrogate for longitude)

   mOutTable[iRowCtr,2:3] <- rEdgeRegionMosaic@data@values[rColVecMosaic[3:4]]  # third column pair:

   mOutTable[iRowCtr,4:5] <- rEdgeRegionCDEM@data@values[rColVecCDEM[3:4]]   # entirely bottom (SRTM part of mosaic) image 

   iRowCtr = iRowCtr + 1  

}

# write the table out as CSV file

message("end of loop - hit key to write output table..")

browser()

write.csv(mOutTable,file="tableForMark4000_5_8.csv",row.names=FALSE)

#

}