Project

General

Profile

« Previous | Next » 

Revision 27c37be1

Added by Rick Reeves over 13 years ago

  • ID 27c37be134db4988f16667468b66183a2a981863

Tested Sunday 5/15 - produces table tableForRic2000_5_8.csv

View differences:

terrain/rscripts/makeImagePairTable.r
1 1 
##############################################################################################
2 2 
#
3 
# makeMarkTable.r
3 
# makeImagePairTable.r
4 4 
#
5 5 
# R script generates table of adjacent ASTER / SRTM / Mosaic Border / CDEM pixel values
6 6 
# to be used to generate plots of 'delta elevation' vs 'pixel pair proximity' and 'elevation'
......
10 10 
# May 4, 2011
11 11 
##############################################################################################
12 12 
#
13 
makeMarkTable <- function()
13 
makeImagePairTable <- function()
14 14 
{
15 15 
require(raster)
16 16 
require(rgdal)
17 17 

  		




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


  		




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


  		




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


  		




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


  		




  
  20
  
    
inputMosaicRaster <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdyFinalBL.tif")


  		




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


  		




  22
  22
  
    
#


  		




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


  		




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


  		




  ......




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


  		




  35
  35
  
    

  		




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


  		




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


  		




  
  37
  
    
eTestAreaExtent <- extent(-135.0,-105.0, 59.995,60.005) # Creates a 12 row subimage


  		




  38
  38
  
    

  		




  39
  39
  
    
# Extract a sub image corresponding to the selected extent.


  		




  40
  40
  
    
# Two different alternatives:


  		




  ......




  52
  52
  
    

  		




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


  		




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


  		




  55
  
  
    

  		




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


  		




  
  55
  
    
extent(rEdgeRegionMosaic) = extent(rEdgeRegionCDEM) # raster package author suggests this to resolve slight extent differences


  		




  
  56
  
    
rEdgeRegionDelta <- rEdgeRegionMosaic - rEdgeRegionCDEM  


  		




  57
  57
  
    

  		




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


  		




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


  		




  60
  60
  
    

  		




  61
  
  
    
nColsToGet <- 4000


  		




  
  61
  
    
nColsToGet <- 2000


  		




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


  		




  
  63
  
    

  		




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


  		




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


  		




  65
  66
  
    

  		




  ......




  79
  80
  
    
{


  		




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


  		




  81
  82
  
    
#browser()


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  
  83
  
    

  		




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


  		




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


  		




  90
  86
  
    
  


  		




  ......




  119
  115
  
    

  		




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


  		




  121
  117
  
    
browser()


  		




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


  		




  
  118
  
    
write.csv(mOutTable,file="/data/project/organisms/rcr/tableForRic2000_5_8.csv",row.names=FALSE)


  		




  123
  119
  
    
#


  		




  124
  120
  
    
}


  		












Also available in:  Unified diff