Project

General

Profile

Download (8.83 KB) Statistics
| Branch: | Revision:
1
##############################################################################################
2
#
3
# SampleDemDiffCols
4
# 
5
# R script generates (or reads in) the CDEM / ASTER-SRTM mosaic 'difference image',
6
# and for randomly-selected one column / two row subimages, computes and saves tbe 
7
# difference between the pixels in the pair to create a distribution of differences.
8
# Three distributions created: North (ASTER CDEM), South i(SRTM/CGIAR), and 
9
# Border (boundary between ASTER and SRTM), Summary statistics are created for each 
10
# distribution.
11
#
12
# Author: Rick Reeves, NCEAS
13
# April 29, 2011
14
##############################################################################################
15
#
16
SampleDemDiffCols <- function()
17
{
18
require(raster)
19
require(rgdal)
20

    
21
#inputRasterMerge <- raster(sFirstImageName)
22
#inputSecondRaster <- raster(sSecondImageName)
23

    
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")
28
#
29
# Difference image for entire merged image takes a while to create, 
30
# so we created it once, now read it back in.
31
#
32
rDeltaWhole <- raster("/data/project/organisms/rcr/ValidateBoundary/DeltaEntireImage.tif")
33
rDeltaWhole@data@values <-getValues(rDeltaWhole)
34

    
35
# Create extent objects used to extract raster subimges. 
36
# The object will be centered along the 60 degree North latitude line,
37
# and have varying depths (number of rows).
38
# The Western Canada study area runs from -135 (west) to -100 (west) longitude,
39
# and 55.0 to 64.00 degrees (north) latitude. 
40
# the ASTER and SRTM/CGIAR image components are merged at the 60 Deg N Latitude line.
41

    
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)
45

    
46
# Extract a sub image corresponding to the selected extent.
47
# Two different alternatives:
48
# The extract() function returns a vector of cell values, 
49
# the crop() function returns a complete raster* object.
50

    
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)
65

    
66
# Important: In order for the image subtraction to work, the extents
67
#            of the two images must be IDENTICAL. I used ArcMap GIS Raster Crop By Mask
68
#            to create subimages with identical extents. 
69

    
70
# Compute the difference image  for the entire study area, and for the region along
71
# the boundary (narrow, maybe 10 pixels either side)
72

    
73
rDeltaEdge <- rEdgeRegionBorder - rEdgeRegionBorderDelta
74

    
75
# Create this image one time, read it in thereafter.
76
#rDeltaWhole <- inputRasterMerge - inputRasterCDEM
77
#writeRaster(rDeltaWhole,filename="DeltaEntireImage.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
78

    
79
# Using the large difference image, compute subimagee statistics for areas
80
# North (ASTER) and South (SRTM) of the boundary. These give us an idea 
81
# re: differences between ASTER and CDEM and CGIAR/SRT and CDEM
82
# what is raster package way of using subscripts to extract? 
83

    
84
# Now, the interesting part: using the boundary difference image, randomly select 
85
# one-degree N-S strips throughout the image, and compare adjacent pixel pairs 
86
# above and below the boundary with pixel pairs straddling the boundary. Subtract  
87
# the pairs, save the collection of (absolute value) of the differences in a vector,
88
# so that we have a population of differences above, below, and straddling the boundary 
89
# line. Compare the populations.
90

    
91
# get a vector of random column index numbers, constrained by column dimension of image
92
# Loop three times, sampling pixel pairs from above, below, across the border
93

    
94
nColsToGet <-20000
95
iDiffVecNorth <- vector(mode="integer",length=nColsToGet)
96
iDiffVecBorder <- vector(mode="integer",length=nColsToGet)
97
iDiffVecSouth <- vector(mode="integer",length=nColsToGet)
98

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

    
101
# Note: initially, sample the same columns in all regions to get a profile.
102
#       other 'sample()' calls can be commented out to sample differenct
103
#       coluns in each 'region'.
104

    
105
# iDiffVecxxxx is a population of differences between adjacent cell pairs. 
106
# Compute iDiffVecNorth/Border/South on either side of border, and across it. 
107
# note that North and South samples taken from larger difference image for 
108
# entire mosaic (sub) image; iDiffBorder taken from the edge region extracted
109
# from the center of the lerger image.
110

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

    
113
colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
114
message("North Sample")
115
#browser()
116
# debug
117
#nColsToGet <- 2
118
#colsToGet <- c(20,100)
119
iFirstRow <- 300
120
iCtr = 1
121
for (iNextCol in colsToGet)
122
{
123
  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
124
  neighborCells <- rDeltaWhole@data@values[rColVec]
125
  iDiffVecNorth[iCtr] <- neighborCells[2] - neighborCells[1]
126
  iCtr = iCtr + 1
127
}
128
#
129
message("Border Sample - different columns")
130
#browser()
131
colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
132
iFirstRow <- 6 # straddle the border of 12 row center section 
133
iCtr = 1
134
for (iNextCol in colsToGet)
135
{
136
  rColVec <- cellFromRowCol(rDeltaEdge,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
137
  neighborCells <- rDeltaEdge@data@values[rColVec]
138
  iDiffVecBorder[iCtr] <- neighborCells[2] - neighborCells[1]
139
  iCtr = iCtr + 1 
140
}
141
#
142
message("South Sample - different columns")
143
#browser()
144
colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
145
iFirstRow <- 3600
146
iCtr = 1
147
for (iNextCol in colsToGet)
148
{
149
  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
150
  neighborCells <- rDeltaWhole@data@values[rColVec]
151
  iDiffVecSouth[iCtr] <- neighborCells[2] - neighborCells[1]
152
  iCtr = iCtr + 1 
153
}
154
# Compute iDiffVecs on either side of border, and across it. 
155
message("Check the cell difference vectors...")
156
#browser()
157

    
158
# summary stats for each population
159

    
160
sNorthSum <- sprintf("ASTER sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
161
                     min(iDiffVecNorth,na.rm=TRUE),median(iDiffVecNorth,na.rm=TRUE),mean(iDiffVecNorth,na.rm=TRUE),
162
                     max(iDiffVecNorth,na.rm=TRUE),var(iDiffVecNorth,na.rm=TRUE),sd(iDiffVecNorth,na.rm=TRUE))
163

    
164
sBorderSum <- sprintf("Border sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
165
                     min(iDiffVecBorder,na.rm=TRUE),median(iDiffVecBorder,na.rm=TRUE),mean(iDiffVecBorder,na.rm=TRUE),
166
                     max(iDiffVecBorder,na.rm=TRUE),var(iDiffVecBorder,na.rm=TRUE),sd(iDiffVecBorder,na.rm=TRUE))
167

    
168
sSouthSum <- sprintf("STRM sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
169
                     min(iDiffVecSouth,na.rm=TRUE),median(iDiffVecSouth,na.rm=TRUE),mean(iDiffVecSouth,na.rm=TRUE),
170
                     max(iDiffVecSouth,na.rm=TRUE),var(iDiffVecSouth,na.rm=TRUE),sd(iDiffVecSouth,na.rm=TRUE))
171
#
172
message(sprintf("statistics for %d N/S adjacent pixel pairs from three mosaic image regions:",nColsToGet))
173
message(sNorthSum)
174
message(sBorderSum)
175
message(sSouthSum)
176

    
177
message("hit key to write output images...")
178
browser()
179

    
180
# Write the extracted subimage and its difference image to disk
181
# For now, use 'gdalinfo' to check image statistics
182

    
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)
186
}
(9-9/19)