Project

General

Profile

Download (10.1 KB) Statistics
| Branch: | Revision:
1
##############################################################################################
2
#
3
# extractDemSubimages
4
# 
5
# R script extracts and saves subimages from various inputs to the 'image boundary analysis';
6
# also produces 'difference images' by subtracting CDEM components from other DEM components.
7
# Output files written as GeoTiff files. 
8
#
9
# Author: Rick Reeves, NCEAS
10
# May 9, 2011
11
##############################################################################################
12
#
13
extractDemSubimages <- function()
14
{
15
require(raster)
16
require(rgdal)
17

    
18
#inputRasterMerge <- raster(sFirstImageName)
19
#inputSecondRaster <- raster(sSecondImageName)
20

    
21
inputRasterAster <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdyASTER_BLEven.tif")
22
inputRasterSRTM <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdySRTM_BLEven.tif")
23
inputRasterMerge <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif")
24
inputRasterCDEM <- raster("/data/project/organisms/rcr/ValidateBoundary/CDEMMosCgiarAsterBdy_BLEven.tif")
25

    
26
inputCdemAster24Row <- raster("/data/project/organisms/rcr/ValidateBoundary/EdgeAsterCDEM24row.tif")
27
inputCdemSrtm24Row <- raster("/data/project/organisms/rcr/ValidateBoundary/EdgeSrtmCDEM24row.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, 60.00009,60.010) # 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, 59.990,60.00) # Creates a 12 row subimage South of border (All SRTM)
45

    
46
eTestAreaExtentAster <- extent(-135.0,-105.0, 60.00009,60.020) # Creates 24 row subimage North of border (All ASTER)
47
eTestAreaExtentBorder <- extent(-135.0,-105.0, 59.99,60.010) # Creates a 24 row subimage centered on border
48
eTestAreaExtentSRTM <- extent(-135.0,-105.0, 59.98,60.00) # Creates a 24 row subimage South of border (All SRTM)
49

    
50

    
51
# Extract a sub image corresponding to the selected extent.
52
# Two different alternatives:
53
# The extract() function returns a vector of cell values, 
54
# the crop() function returns a complete raster* object.
55
message("extracts")
56
browser()
57
vEdgeRegionAster <- extract(inputRasterAster,eTestAreaExtentAster)
58
rEdgeRegionAster <- crop(inputRasterAster,eTestAreaExtentAster)
59
vEdgeRegionAsterDelta <- extract(rDeltaWhole,eTestAreaExtentAster)
60
rEdgeRegionAsterDelta <- crop(rDeltaWhole,eTestAreaExtentAster)
61
#vEdgeRegionAsterCDEM <- extract(inputRasterCDEM,eTestAreaExtentAster)
62
#rEdgeRegionAsterCDEM <- crop(inputRasterCDEM,eTestAreaExtentAster)
63

    
64
vEdgeRegionBorder <- extract(inputRasterMerge,eTestAreaExtentBorder)
65
rEdgeRegionBorder <- crop(inputRasterMerge,eTestAreaExtentBorder)
66
vEdgeRegionBorderDelta <- extract(rDeltaWhole,eTestAreaExtentBorder)
67
rEdgeRegionBorderDelta <- crop(rDeltaWhole,eTestAreaExtentBorder)
68

    
69
vEdgeRegionSRTM <- extract(inputRasterSRTM,eTestAreaExtentSRTM)
70
rEdgeRegionSRTM <- crop(inputRasterSRTM,eTestAreaExtentSRTM)
71
vEdgeRegionSRTMDelta <- extract(rDeltaWhole,eTestAreaExtentSRTM)
72
rEdgeRegionSRTMDelta <- crop(rDeltaWhole,eTestAreaExtentSRTM)
73
#vEdgeRegionSrtmCDEM <- extract(inputRasterCDEM,eTestAreaExtentSRTM)
74
#rEdgeRegionSrtmCDEM <- crop(inputRasterCDEM,eTestAreaExtentSRTM)
75
message("differences")
76
browser()
77
# 
78
# write these files to disk, get stats.
79
#
80
rDeltaEdgeAster <- rEdgeRegionAster - inputCdemAster24Row
81
rDeltaEdgeBorder <- rEdgeRegionBorder - rEdgeRegionBorderDelta
82
rDeltaEdgeSRTM <- rEdgeRegionSRTM - rEdgeRegionSrtmCDEM
83

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

    
88
# Computeexter the difference image  for the entire study area, and for the region along
89
# the boundary (narrow, maybe 10 pixels either side)
90

    
91
# rDeltaEdge <- rEdgeRegionFirst - rEdgeRegionSecond
92
# May 10 2011 note: One of these statements generates the 'different origin' fatal error. 
93
# my workaround was to create ALL of the subimages used here with gdalwarp, and then 
94
# read them in and subtract them. 
95
# see other R script: OnlyImage Differences.r
96
# I will create a reproducable example, send to Raster authors. 
97

    
98
rDeltaEdgeAster <- rEdgeRegionAster - inputCdemAster24Row
99
rDeltaEdgeSRTM <- rEdgeRegionSRTM - inputCdemSrtm24Row
100

    
101
#rDeltaWhole <- inputRasterMerge - inputRasterCDEM
102
#writeRaster(rDeltaWhole,filename="DeltaEntireImage.tif",format="GTiff",datatype="INT2S",overwrite=TRUE)
103

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

    
109
# Now, the interesting part: using the boundary difference image, randomly select 
110
# one-degree N-S strips throughout the image, and compare adjacent pixel pairs 
111
# above and below the boundary with pixel pairs straddling the boundary. Subtract  
112
# the pairs, save the collection of (absolute value) of the differences in a vector,
113
# so that we have a population of differences above, below, and straddling the boundary 
114
# line. Compare the populations.
115

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

    
119
nColsToGet <-20000
120
iDiffVecNorth <- vector(mode="integer",length=nColsToGet)
121
iDiffVecBorder <- vector(mode="integer",length=nColsToGet)
122
iDiffVecSouth <- vector(mode="integer",length=nColsToGet)
123

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

    
126
# Note: initially, sample the same columns in all regions to get a profile.
127
#       other 'sample()' calls can be commented out to sample differenct
128
#       coluns in each 'region'.
129

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

    
136
# Remember, we are sampling a PAIR of pixels (same column from two adjacent rows)
137
# NOTR 5/9: need to rethink these sampling loops for this new version. Comment for now!
138
#colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
139
#message("North Sample")
140
#browser()
141
# debug
142
#nColsToGet <- 2
143
#colsToGet <- c(20,100)
144
#iFirstRow <- 300
145
#iCtr = 1
146
#for (iNextCol in colsToGet)
147
#{
148
#  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
149
#  neighborCells <- rDeltaWhole@data@values[rColVec]
150
#  iDiffVecNorth[iCtr] <- neighborCells[2] - neighborCells[1]
151
#  iCtr = iCtr + 1
152
#}
153
#
154
#message("Border Sample - different columns")
155
#browser()
156
#colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
157
#iFirstRow <- 6 # straddle the border of 12 row center section 
158
#iCtr = 1
159
#for (iNextCol in colsToGet)
160
#{
161
#  rColVec <- cellFromRowCol(rDeltaEdge,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
162
#  neighborCells <- rDeltaEdge@data@values[rColVec]
163
#  iDiffVecBorder[iCtr] <- neighborCells[2] - neighborCells[1]
164
#  iCtr = iCtr + 1 
165
#}
166
##
167
#message("South Sample - different columns")
168
##browser()
169
#colsToGet <-sample(1:inputRasterMerge@ncols,nColsToGet)
170
#iFirstRow <- 3600
171
#iCtr = 1
172
#for (iNextCol in colsToGet)
173
#{
174
#  rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol)
175
#  neighborCells <- rDeltaWhole@data@values[rColVec]
176
#  iDiffVecSouth[iCtr] <- neighborCells[2] - neighborCells[1]
177
#  iCtr = iCtr + 1 
178
#}
179
## Compute iDiffVecs on either side of border, and across it. 
180
#message("Check the cell difference vectors...")
181
##browser()
182

    
183
# summary stats for each population
184

    
185
#sNorthSum <- sprintf("ASTER sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
186
#                     min(iDiffVecNorth,na.rm=TRUE),median(iDiffVecNorth,na.rm=TRUE),mean(iDiffVecNorth,na.rm=TRUE),
187
#                     max(iDiffVecNorth,na.rm=TRUE),var(iDiffVecNorth,na.rm=TRUE),sd(iDiffVecNorth,na.rm=TRUE))
188
#
189
#sBorderSum <- sprintf("Border sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
190
#                    min(iDiffVecBorder,na.rm=TRUE),median(iDiffVecBorder,na.rm=TRUE),mean(iDiffVecBorder,na.rm=TRUE),
191
#                    max(iDiffVecBorder,na.rm=TRUE),var(iDiffVecBorder,na.rm=TRUE),sd(iDiffVecBorder,na.rm=TRUE))
192

    
193
#sSouthSum <- sprintf("STRM sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
194
#                     min(iDiffVecSouth,na.rm=TRUE),median(iDiffVecSouth,na.rm=TRUE),mean(iDiffVecSouth,na.rm=TRUE),
195
#                     max(iDiffVecSouth,na.rm=TRUE),var(iDiffVecSouth,na.rm=TRUE),sd(iDiffVecSouth,na.rm=TRUE))
196
##
197
#message(sprintf("statistics for %d N/S adjacent pixel pairs from three mosaic image regions:",nColsToGet))
198
#message(sNorthSum)
199
#message(sBorderSum)
200
#message(sSouthSum)
201

    
202
message("hit key to write output images...")
203
browser()
204

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

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

    
210
}
(10-10/19)