Revision 5c833657
Added by Rick Reeves over 13 years ago
- ID 5c833657f2f5a7b0e8678d380ff2b51940fa7e6d
| terrain/rscripts/SampleDemDiffCols.r | ||
|---|---|---|
| 13 | 13 |
# April 29, 2011 |
| 14 | 14 |
############################################################################################## |
| 15 | 15 |
# |
| 16 |
#ExtractDemRasterSubimages <- function(sFirstImageName,sSecondImageName) |
|
| 17 | 16 |
SampleDemDiffCols <- function() |
| 18 | 17 |
{
|
| 19 | 18 |
require(raster) |
| 20 | 19 |
require(rgdal) |
| 21 | 20 |
|
| 22 |
#inputFirstRaster <- raster(sFirstImageName)
|
|
| 21 |
#inputRasterMerge <- raster(sFirstImageName)
|
|
| 23 | 22 |
#inputSecondRaster <- raster(sSecondImageName) |
| 24 | 23 |
|
| 25 |
inputFirstRaster <- raster("mergeCgiarAsterBdyTuesdayClip.tif")
|
|
| 26 |
inputSecondRaster <- raster("CDemMosTuesdayClipMergeSpace.tif")
|
|
| 24 |
inputRasterAster <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")
|
|
| 25 |
inputRasterSRTM <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")
|
|
| 26 |
inputRasterMerge <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")
|
|
| 27 |
inputRasterCDEM <- raster("/data/project/organisms/rcr/ValidateBoundary/CDemMosTuesdayClipMergeSpace.tif")
|
|
| 27 | 28 |
# |
| 28 | 29 |
# Difference image for entire merged image takes a while to create, |
| 29 | 30 |
# so we created it once, now read it back in. |
| 30 | 31 |
# |
| 31 |
rDeltaWhole <- raster("DeltaEntireImage.tif")
|
|
| 32 |
rDeltaWhole <- raster("/data/project/organisms/rcr/ValidateBoundary/DeltaEntireImage.tif")
|
|
| 32 | 33 |
rDeltaWhole@data@values <-getValues(rDeltaWhole) |
| 33 | 34 |
|
| 34 | 35 |
# Create extent objects used to extract raster subimges. |
| ... | ... | |
| 38 | 39 |
# and 55.0 to 64.00 degrees (north) latitude. |
| 39 | 40 |
# the ASTER and SRTM/CGIAR image components are merged at the 60 Deg N Latitude line. |
| 40 | 41 |
|
| 41 |
#eTestAreaExtent <- extent(-135.2,-100.2, 59.997,60.001) # Creates a 5 row subimage |
|
| 42 |
eTestAreaExtent <- extent(-135.2,-100.2, 59.995,60.005) # Creates a 12 row subimage |
|
| 42 |
eTestAreaExtentAster <- extent(-135.0,-105.0, 59.990,60.00) # Creates 12 row subimage North of border (All ASTER) |
|
| 43 |
eTestAreaExtentBorder <- extent(-135.0,-105.0, 59.995,60.005) # Creates a 12 row subimage centered on border |
|
| 44 |
eTestAreaExtentSRTM <- extent(-135.0,-105.0, 60.00,60.010) # Creates a 12 row subimage South of border (All SRTM) |
|
| 43 | 45 |
|
| 44 | 46 |
# Extract a sub image corresponding to the selected extent. |
| 45 | 47 |
# Two different alternatives: |
| 46 | 48 |
# The extract() function returns a vector of cell values, |
| 47 | 49 |
# the crop() function returns a complete raster* object. |
| 48 | 50 |
|
| 49 |
vEdgeRegionFirst <- extract(inputFirstRaster,eTestAreaExtent) |
|
| 50 |
rEdgeRegionFirst <- crop(inputFirstRaster,eTestAreaExtent) |
|
| 51 |
vEdgeRegionSecond <- extract(inputSecondRaster,eTestAreaExtent) |
|
| 52 |
rEdgeRegionSecond <- crop(inputSecondRaster,eTestAreaExtent) |
|
| 51 |
vEdgeRegionAster <- extract(inputRasterMerge,eTestAreaExtentAster) |
|
| 52 |
rEdgeRegionAster <- crop(inputRasterMerge,eTestAreaExtentAster) |
|
| 53 |
vEdgeRegionAsterDelta <- extract(rDeltaWhole,eTestAreaExtentAster) |
|
| 54 |
rEdgeRegionAsterDelta <- crop(rDeltaWhole,eTestAreaExtentAster) |
|
| 55 |
|
|
| 56 |
vEdgeRegionBorder <- extract(inputRasterCDEM,eTestAreaExtentBorder) |
|
| 57 |
rEdgeRegionBorder <- crop(inputRasterMerge,eTestAreaExtentBorder) |
|
| 58 |
vEdgeRegionBorderDelta <- extract(rDeltaWhole,eTestAreaExtentBorder) |
|
| 59 |
rEdgeRegionBorderDelta <- crop(rDeltaWhole,eTestAreaExtentBorder) |
|
| 60 |
|
|
| 61 |
vEdgeRegionSrtm <- extract(inputRasterMerge,eTestAreaExtentSRTM) |
|
| 62 |
rEdgeRegionSrtm <- crop(inputRasterMerge,eTestAreaExtentSRTM) |
|
| 63 |
vEdgeRegionSRTMDelta <- extract(rDeltaWhole,eTestAreaExtentSRTM) |
|
| 64 |
rEdgeRegionSRTMDelta <- crop(rDeltaWhole,eTestAreaExtentSRTM) |
|
| 53 | 65 |
|
| 54 | 66 |
# Important: In order for the image subtraction to work, the extents |
| 55 | 67 |
# of the two images must be IDENTICAL. I used ArcMap GIS Raster Crop By Mask |
| ... | ... | |
| 58 | 70 |
# Compute the difference image for the entire study area, and for the region along |
| 59 | 71 |
# the boundary (narrow, maybe 10 pixels either side) |
| 60 | 72 |
|
| 61 |
rDeltaEdge <- rEdgeRegionFirst - rEdgeRegionSecond
|
|
| 73 |
rDeltaEdge <- rEdgeRegionBorder - rEdgeRegionBorderDelta
|
|
| 62 | 74 |
|
| 63 | 75 |
# Create this image one time, read it in thereafter. |
| 64 |
#rDeltaWhole <- inputFirstRaster - inputSecondRaster
|
|
| 76 |
#rDeltaWhole <- inputRasterMerge - inputRasterCDEM
|
|
| 65 | 77 |
#writeRaster(rDeltaWhole,filename="DeltaEntireImage.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
| 66 | 78 |
|
| 67 | 79 |
# Using the large difference image, compute subimagee statistics for areas |
| ... | ... | |
| 79 | 91 |
# get a vector of random column index numbers, constrained by column dimension of image |
| 80 | 92 |
# Loop three times, sampling pixel pairs from above, below, across the border |
| 81 | 93 |
|
| 82 |
nColsToGet <- 1000
|
|
| 94 |
nColsToGet <-20000
|
|
| 83 | 95 |
iDiffVecNorth <- vector(mode="integer",length=nColsToGet) |
| 84 | 96 |
iDiffVecBorder <- vector(mode="integer",length=nColsToGet) |
| 85 | 97 |
iDiffVecSouth <- vector(mode="integer",length=nColsToGet) |
| ... | ... | |
| 98 | 110 |
|
| 99 | 111 |
# Remember, we are sampling a PAIR of pixels (same column from two adjacent rows) |
| 100 | 112 |
|
| 101 |
colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet)
|
|
| 113 |
colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
|
|
| 102 | 114 |
message("North Sample")
|
| 103 | 115 |
#browser() |
| 104 | 116 |
# debug |
| ... | ... | |
| 110 | 122 |
{
|
| 111 | 123 |
rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol) |
| 112 | 124 |
neighborCells <- rDeltaWhole@data@values[rColVec] |
| 113 |
iDiffVecNorth[iCtr] <- neighborCells[1] - neighborCells[2]
|
|
| 125 |
iDiffVecNorth[iCtr] <- neighborCells[2] - neighborCells[1]
|
|
| 114 | 126 |
iCtr = iCtr + 1 |
| 115 | 127 |
} |
| 116 | 128 |
# |
| 117 |
message("Border Sample")
|
|
| 129 |
message("Border Sample - different columns")
|
|
| 118 | 130 |
#browser() |
| 119 |
#colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet)
|
|
| 131 |
colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
|
|
| 120 | 132 |
iFirstRow <- 6 # straddle the border of 12 row center section |
| 121 | 133 |
iCtr = 1 |
| 122 | 134 |
for (iNextCol in colsToGet) |
| 123 | 135 |
{
|
| 124 | 136 |
rColVec <- cellFromRowCol(rDeltaEdge,iFirstRow:(iFirstRow+1),iNextCol:iNextCol) |
| 125 | 137 |
neighborCells <- rDeltaEdge@data@values[rColVec] |
| 126 |
iDiffVecBorder[iCtr] <- neighborCells[1] - neighborCells[2]
|
|
| 138 |
iDiffVecBorder[iCtr] <- neighborCells[2] - neighborCells[1]
|
|
| 127 | 139 |
iCtr = iCtr + 1 |
| 128 | 140 |
} |
| 129 | 141 |
# |
| 130 |
message("South Sample")
|
|
| 142 |
message("South Sample - different columns")
|
|
| 131 | 143 |
#browser() |
| 132 |
#colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet)
|
|
| 144 |
colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
|
|
| 133 | 145 |
iFirstRow <- 3600 |
| 134 | 146 |
iCtr = 1 |
| 135 | 147 |
for (iNextCol in colsToGet) |
| 136 | 148 |
{
|
| 137 | 149 |
rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol) |
| 138 | 150 |
neighborCells <- rDeltaWhole@data@values[rColVec] |
| 139 |
iDiffVecSouth[iCtr] <- neighborCells[1] - neighborCells[2]
|
|
| 151 |
iDiffVecSouth[iCtr] <- neighborCells[2] - neighborCells[1]
|
|
| 140 | 152 |
iCtr = iCtr + 1 |
| 141 | 153 |
} |
| 142 | 154 |
# Compute iDiffVecs on either side of border, and across it. |
| ... | ... | |
| 164 | 176 |
|
| 165 | 177 |
message("hit key to write output images...")
|
| 166 | 178 |
browser() |
| 179 |
|
|
| 167 | 180 |
# Write the extracted subimage and its difference image to disk |
| 168 | 181 |
# For now, use 'gdalinfo' to check image statistics |
| 169 | 182 |
|
| 170 |
writeRaster(rEdgeRegionFirst,filename="EdgeRegionFirst.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
| 171 |
|
|
| 172 |
writeRaster(rEdgeRegionAsterCgiar,filename="EdgeRegionSecond.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
| 173 |
|
|
| 174 |
writeRaster(rDelta,filename="EdgeRegionDelta.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
| 183 |
writeRaster(rEdgeRegionBorder,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionFirstDC.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
| 184 |
writeRaster(rEdgeRegionBorderDelta,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionSecondDC.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
| 185 |
writeRaster(rDeltaEdge,filename="/data/project/organisms/rcr/ValidateBoundary/EdgeRegionDeltaDC.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
| 175 | 186 |
} |
Also available in: Unified diff
Version runs correctly, though the CDEM edge input image might not be the latest and best. Recheck inputs