Revision ad35c236
Added by Rick Reeves over 13 years ago
- ID ad35c236f35d04cde989aaf17ee6e683a533a213
terrain/rscripts/makeMarkTable.r | ||
---|---|---|
15 | 15 |
require(raster) |
16 | 16 |
require(rgdal) |
17 | 17 |
|
18 |
#inputFirstRaster <- raster(sFirstImageName) |
|
19 |
#inputSecondRaster <- raster(sSecondImageName) |
|
20 |
|
|
21 | 18 |
inputCgiarRaster <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdySRTM_BL.tif") |
22 | 19 |
inputAsterRaster <- raster("/data/project/organisms/rcr/AsterCgiarMerge/mergeCgiarAsterBdyASTER_BL.tif") |
23 | 20 |
inputMosaicRaster <- raster("/data/project/organisms/rcr/ValidateBoundary/mergeCgiarAsterBdyTuesdayClip.tif") |
... | ... | |
44 | 41 |
# The extract() function returns a vector of cell values, |
45 | 42 |
# the crop() function returns a complete raster* object. |
46 | 43 |
|
47 |
#vEdgeRegionMosaic <- extract(inputMosaicRaster,eTestAreaExtent) |
|
48 |
#vEdgeRegionCDEM <- extract(inputCDEMRaster,eTestAreaExtent) |
|
49 |
|
|
50 | 44 |
rEdgeRegionAster <- crop(inputAsterRaster,eTestAreaExtent) |
51 | 45 |
rEdgeRegionCgiar <- crop(inputCgiarRaster,eTestAreaExtent) |
52 | 46 |
rEdgeRegionMosaic <- crop(inputMosaicRaster,eTestAreaExtent) |
... | ... | |
61 | 55 |
|
62 | 56 |
rEdgeRegionDelta <- rEdgeRegionMosaic - rEdgeRegionCDEM # not used in this version (yet) |
63 | 57 |
|
64 |
# Create this image one time, read it in thereafter. |
|
65 |
#rDeltaWhole <- inputMosaicRaster - inputCDEMRaster |
|
66 |
#writeRaster(rDeltaWhole,filename="DeltaMosaicCDEMSubmage.tif",format="GTiff",datatype="INT2S",overwrite=TRUE) |
|
67 |
|
|
68 |
# Using the large difference image, compute subimagee statistics for areas |
|
69 |
# North (ASTER) and South (SRTM) of the boundary. These give us an idea |
|
70 |
# re: differences between ASTER and CDEM and CGIAR/SRT and CDEM |
|
71 |
# what is raster package way of using subscripts to extract? |
|
72 |
|
|
73 |
# Now, the interesting part: using the boundary difference image, randomly select |
|
74 |
# one-degree N-S strips throughout the image, and compare adjacent pixel pairs |
|
75 |
# above and below the boundary with pixel pairs straddling the boundary. Subtract |
|
76 |
# the pairs, save the collection of (absolute value) of the differences in a vector, |
|
77 |
# so that we have a population of differences above, below, and straddling the boundary |
|
78 |
# line. Compare the populations. |
|
79 |
|
|
80 | 58 |
# get a vector of random column index numbers, constrained by column dimension of image |
81 | 59 |
# Loop three times, sampling pixel pairs from above, below, across the border |
82 | 60 |
|
83 |
nColsToGet <- 1000
|
|
61 |
nColsToGet <- 4000
|
|
84 | 62 |
iDiffVecNorth <- vector(mode="integer",length=nColsToGet) |
85 | 63 |
iDiffVecBorder <- vector(mode="integer",length=nColsToGet) |
86 | 64 |
iDiffVecSouth <- vector(mode="integer",length=nColsToGet) |
87 | 65 |
|
88 |
#colsToGet <-sample(1:50,nColsToGet) |
|
89 |
|
|
90 |
# Note: initially, sample the same columns in all regions to get a profile. |
|
91 |
# other 'sample()' calls can be commented out to sample differenct |
|
92 |
# coluns in each 'region'. |
|
93 |
|
|
94 |
# iDiffVecxxxx is a population of differences between adjacent cell pairs. |
|
95 |
# Compute iDiffVecNorth/Border/South on either side of border, and across it. |
|
96 |
# note that North and South samples taken from larger difference image for |
|
97 |
# entire mosaic (sub) image; iDiffBorder taken from the edge region extracted |
|
98 |
# from the center of the lerger image. |
|
99 |
|
|
100 |
# Remember, we are sampling a PAIR of pixels (same column from two adjacent rows) |
|
101 |
|
|
102 |
#browser() |
|
103 |
# debug |
|
104 |
#nColsToGet <- 2 |
|
105 |
#colsToGet <- c(20,100) |
|
106 |
#iFirstRow <- 300 |
|
107 |
#iCtr = 1 |
|
108 |
#for (iNextCol in colsToGet) |
|
109 |
#{ |
|
110 |
# rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol) |
|
111 |
# neighborCells <- rDeltaWhole@data@values[rColVec] |
|
112 |
# iDiffVecNorth[iCtr] <- neighborCells[2] - neighborCells[1] |
|
113 |
#iCtr = iCtr + 1 |
|
114 |
#} |
|
115 |
# |
|
116 |
|
|
117 | 66 |
# Here is the output matrix, (nColSamples * 3) rows, four columns |
118 | 67 |
colsToGet <-sample(1:inputMosaicRaster@ncols,nColsToGet) |
119 |
message("North Sample") |
|
120 | 68 |
|
121 | 69 |
mOutTable <- matrix(nrow=(nColsToGet * 3), ncol=5) |
122 | 70 |
colnames(mOutTable) <- c("ColumnID","elevNorth","elevSouth","cdemNorth","cdemSouth") |
123 | 71 |
# in this difference image, the border edge occurs at column 6 |
124 |
message("Border Sample - different columns") |
|
125 |
#browser() |
|
126 |
#colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet) |
|
127 |
iFirstRow <- 4 # two rows before the border |
|
128 |
iLastRow <- 7 # two rows after the border |
|
72 |
|
|
73 |
# NOTE Wed eve: A better start and end for the border: firstRow=5, lastRow = 8.... |
|
74 |
iFirstRow <- 5 #4 # two rows before the border |
|
75 |
iLastRow <- 8 #7 # two rows after the border |
|
129 | 76 |
iRowCtr = 1 # points to latest output table being written |
77 |
message(sprintf("starting col sample loop")) |
|
130 | 78 |
for (iNextCol in colsToGet) |
131 | 79 |
{ |
132 |
rColVecAster <- cellFromRowCol(rEdgeRegionAster,iFirstRow:(iLastRow),iNextCol:iNextCol) |
|
133 |
rColVecCgiar <- cellFromRowCol(rEdgeRegionCgiar,iFirstRow:(iLastRow),iNextCol:iNextCol) |
|
134 |
rColVecMosaic <- cellFromRowCol(rEdgeRegionMosaic,iFirstRow:(iLastRow),iNextCol:iNextCol) |
|
135 |
rColVecCDEM <- cellFromRowCol(rEdgeRegionCDEM,iFirstRow:(iLastRow),iNextCol:iNextCol) |
|
80 |
#message(sprintf("in col sample loop getting col sample %d - hit key..",iNextCol)) |
|
81 |
#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) |
|
88 |
rColVecMosaic <- cellFromRowCol(rEdgeRegionMosaic,iFirstRow:(iLastRow),iNextCol:iNextCol) |
|
89 |
rColVecCDEM <- cellFromRowCol(rEdgeRegionCDEM,iFirstRow:(iLastRow),iNextCol:iNextCol) |
|
136 | 90 |
|
137 | 91 |
# Split the column vector into the pairs that Mark requested |
138 | 92 |
# For each column sampled, the output table has three rows: |
139 | 93 |
# |
140 | 94 |
# Row 1: from (input) Aster layer: Two pixels above border |
141 |
# Row 2: from (input) Mosaic layer: Two pixels straddling border |
|
95 |
# Row 2: from (input) Mosaic layer: Two pixels straddling border rDeltaWhole@data@values[rColVec]
|
|
142 | 96 |
# Row 3: from (input) Cgiar layer: Two pixels below border |
143 | 97 |
# |
144 | 98 |
# each with four columns, arranged into two column pairs: |
145 | 99 |
# |
146 | 100 |
# North Pixel and South Pixel elevation, North Pixel and South Pixel CDEM (baseline) |
147 | 101 |
# |
148 |
mOutTable[iRowCtr][1] <- iNextCol # The (radomly sampled) col ID (surrogate for longitude) |
|
149 |
mOutTable[iRowCtr][2:3] <- rColVecAster[1:2] # First column pair from extracted vector: |
|
150 |
mOutTable[iRowCtr][4:5] <- rColVecCDEM[1:2] # entirely top (ASTER) image |
|
151 |
iRowCtr = iRowCtr + 1 |
|
152 |
# |
|
153 |
mOutTable[iRowCtr][1] <- iNextCol # The (radomly sampled) col ID (surrogate for longitude) |
|
154 |
mOutTable[iRowCtr][2:3] <- rColVecMosaic[3:4] # Second column pair: |
|
155 |
mOutTable[iRowCtr][4:5] <- rColVecCDEM[3:4] # straddles border region (from Mosaic image) |
|
156 |
iRowCtr = iRowCtr + 1 |
|
102 |
mOutTable[iRowCtr,1] <- iNextCol # The (radomly sampled) col ID (surrogate for longitude) |
|
103 |
mOutTable[iRowCtr,2:3] <- rEdgeRegionMosaic@data@values[rColVecMosaic[1:2]] # First column pair from extracted vector: |
|
104 |
mOutTable[iRowCtr,4:5] <- rEdgeRegionCDEM@data@values[rColVecCDEM[1:2]] # entirely top (ASTER part of mosaic) image |
|
105 |
iRowCtr = iRowCtr + 1 |
|
157 | 106 |
# |
158 |
mOutTable[iRowCtr][1] <- iNextCol # The (radomly sampled) col ID (surrogate for longitude) |
|
159 |
mOutTable[iRowCtr][2:3] <- rColVecCgiar[5:6] # third column pair: |
|
160 |
mOutTable[iRowCtr][4:5] <- rColVecCDEM[5:6] # entirely bottom (SRTM) image |
|
161 |
iRowCtr = iRowCtr + 1 |
|
162 |
# |
|
107 |
mOutTable[iRowCtr,1] <- iNextCol # The (radomly sampled) col ID (surrogate for longitude) |
|
108 |
mOutTable[iRowCtr,2:3] <- rEdgeRegionMosaic@data@values[rColVecMosaic[2:3]] # Second column pair: |
|
109 |
mOutTable[iRowCtr,4:5] <- rEdgeRegionCDEM@data@values[rColVecCDEM[2:3]] # straddles border region |
|
110 |
iRowCtr = iRowCtr + 1 |
|
111 |
# |
|
112 |
mOutTable[iRowCtr,1] <- iNextCol # The (radomly sampled) col ID (surrogate for longitude) |
|
113 |
mOutTable[iRowCtr,2:3] <- rEdgeRegionMosaic@data@values[rColVecMosaic[3:4]] # third column pair: |
|
114 |
mOutTable[iRowCtr,4:5] <- rEdgeRegionCDEM@data@values[rColVecCDEM[3:4]] # entirely bottom (SRTM part of mosaic) image |
|
115 |
iRowCtr = iRowCtr + 1 |
|
163 | 116 |
} |
164 |
# |
|
165 |
# write the table out as CSV file |
|
166 |
# TODO: add the column ID (proximity indicator |
|
167 |
message("hit key to write output table..") |
|
168 |
browser() |
|
169 |
writeCSV(mOutTable,file="tableForMark.csv",row.names=FALSE) |
|
170 |
# |
|
171 |
#message("South Sample - different columns") |
|
172 |
#browser() |
|
173 |
#colsToGet <-sample(1:inputFirstRaster@ncols,nColsToGet) |
|
174 |
#iFirstRow <- 3600 |
|
175 |
#iCtr = 1 |
|
176 |
#for (iNextCol in colsToGet) |
|
177 |
#{ |
|
178 |
# rColVec <- cellFromRowCol(rDeltaWhole,iFirstRow:(iFirstRow+1),iNextCol:iNextCol) |
|
179 |
# neighborCells <- rDeltaWhole@data@values[rColVec] |
|
180 |
# iDiffVecSouth[iCtr] <- neighborCells[2] - neighborCells[1] |
|
181 |
# iCtr = iCtr + 1 |
|
182 |
#} |
|
183 |
# Compute iDiffVecs on either side of border, and across it. |
|
184 |
message("Check the cell difference vectors...") |
|
185 |
#browser() |
|
186 |
|
|
187 |
# summary stats for each population |
|
188 |
|
|
189 |
#sNorthSum <- sprintf("ASTER sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f", |
|
190 |
# min(iDiffVecNorth,na.rm=TRUE),median(iDiffVecNorth,na.rm=TRUE),mean(iDiffVecNorth,na.rm=TRUE), |
|
191 |
# max(iDiffVecNorth,na.rm=TRUE),var(iDiffVecNorth,na.rm=TRUE),sd(iDiffVecNorth,na.rm=TRUE)) |
|
192 | 117 |
|
193 |
sBorderSum <- sprintf("Border sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f", |
|
194 |
min(iDiffVecBorder,na.rm=TRUE),median(iDiffVecBorder,na.rm=TRUE),mean(iDiffVecBorder,na.rm=TRUE), |
|
195 |
max(iDiffVecBorder,na.rm=TRUE),var(iDiffVecBorder,na.rm=TRUE),sd(iDiffVecBorder,na.rm=TRUE)) |
|
118 |
# write the table out as CSV file |
|
196 | 119 |
|
197 |
#sSouthSum <- sprintf("STRM sample summary: Min: %f / Median: %d / Mean: %f / Max: %f / Variance: %f sDev: %f",
|
|
198 |
# min(iDiffVecSouth,na.rm=TRUE),median(iDiffVecSouth,na.rm=TRUE),mean(iDiffVecSouth,na.rm=TRUE),
|
|
199 |
# max(iDiffVecSouth,na.rm=TRUE),var(iDiffVecSouth,na.rm=TRUE),sd(iDiffVecSouth,na.rm=TRUE))
|
|
120 |
message("end of loop - hit key to write output table..")
|
|
121 |
browser()
|
|
122 |
write.csv(mOutTable,file="tableForMark4000_5_8.csv",row.names=FALSE)
|
|
200 | 123 |
# |
201 |
message(sprintf("statistics for %d N/S adjacent pixel pairs from three mosaic image regions:",nColsToGet)) |
|
202 |
#message(sNorthSum) |
|
203 |
message(sBorderSum) |
|
204 |
#mssage(sSouthSum) |
|
205 |
|
|
206 | 124 |
} |
Also available in: Unified diff
Producing valid table, but from column samples taken entirely from mosaic. Next step: take them from Aster and Cgiar images.