Project

General

Profile

Download (7.86 KB) Statistics
| Branch: | Revision:
1

    
2
# DEMBdyFixExpDecay - Implements Exponential Decay method to 'blend' DEM boundary edges
3
# 
4
# General form: decayFactor = e(-(d * DKc))
5
#                 where d = distance interval from origin (in this case, an image row index)
6
#                       DKc = empirically derived constant
7
# 
8
#
9
# Inputs: 
10
#   1) Blended Aster/SRTM DEM mosaic image.
11
#   2) 'SRTM-Only' subimage corresponding to the NORTHernmost part (adjacent to the 60 degree Lat 
12
#     boundary) of the SRTM component of the mosaic.
13
#   3) 'Aster-Only' subimage corresponding to the SOUTHernmost part (adjacent to the 60 degree Lat 
14
#     boundary) of the ASTER component of the mosaic.
15
#   4) Other parameters, as they are needed
16
#
17
# To run:
18

    
19
 
20
#  2) > source("DEMBdyFixExpDecay.r")
21
#  3) > DEMBdyFixExpDecay()
22
#
23
# Author: Rick Reeves, NCEAS
24
# June, 2011
25

    
26
##############################################################################################
27
DEMBdyFixExpDecay <- function()
28
{
29
  require(raster)
30
  require(rgdal)
31
  
32
# Read the mosaic image and component sub-images:
33
#   - Aster/SRTM mosaic image (with boundary edge to be fixed)
34
#   - 'Aster-only' subimage adjacent to Southern edge of 60 degree bounls -s dary,
35
#     sampled into a grid with same extent and spatial resolution as image mosaic
36
#   - 'Aster-only' subimage adjacent to Northern edge of 60 degree boundary,
37
#     sampled into a grid with same extent and spatial resolution as image mosaic
38
    
39
     inMosaic      <- raster("./AsterSrtmBoth3ArcSecSub.tif")
40
     
41
# Use these objects to get the extents of the subsets within the mosaic
42

    
43
     inAsterNorth  <- raster("./AsterBdyTestAbove60Sub.tif")     
44
     inAsterSouth  <- raster("./AsterBdyTestBelow60Sub.tif")     
45
     
46
# Better, read them into a rasterStack
47

    
48
     mosaicLayers <- stack("./AsterSRTMBothTestSub.tif",
49
                           "./AsterBdyTestFGAbove60Sub.tif",
50
                           "./AsterBdyTestFGBelow60Sub.tif") 
51
# 
52
# Create copy of input raster that we will use to create a 'fixed boundary' iamge
53
# Even SIMPLER, according to raster documentation
54
#
55
     outMosaic <- raster(mosaicLayers[[1]])
56
 
57
     sFixedRasterFile <- "TestOutputRasterBdyFix.tif"
58
     
59
# First, get the extent of the 'below 60' sub image in the (big) input raster
60

    
61
     northAsterEx <- extent(inAsterNorth)
62
     southAsterEx <- extent(inAsterSouth)
63
     
64
# Get the values from the mosaic for this extent
65

    
66
     southCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],southAsterEx)
67
     northCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],northAsterEx)
68

    
69
# Within the large input raster, we need the index of the first row of the 'Below 60 Aster' subimage.
70
# Our plan: to replace this portion of the input mosaic with a linear combination of the original
71
# input mosaic COPY (outMosaic) and the 'below 60' raster.
72

    
73
     firstNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[1])   
74
     lastNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[(length(northCellsOfInterest) - 1)])
75
     northNrowsToProcess <- nrow(inAsterNorth)
76
  
77
     firstSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[1])   
78
     lastSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[(length(southCellsOfInterest) - 1)])
79
     southNrowsToProcess <- nrow(inAsterSouth)
80

    
81
#  create the output raster by copying the input file with 'block copy'
82
     
83
message("Create output (fixed) mosaic")
84
#browser()
85
     bs <- blockSize(mosaicLayers[[1]]) 
86
     outMosaic <- writeStart(outMosaic,sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
87
     for (iCtr in 1 : bs$n)
88
     {
89
        mosaicVals <- getValues(mosaicLayers[[1]],row=bs$row[iCtr],nrows=bs$nrows[iCtr])
90
        writeValues(outMosaic,mosaicVals,bs$row[iCtr])
91
        message(sprintf(".....processed input mosaic block %d",iCtr))        
92
     }
93
     outMosaic <- writeStop(outMosaic)
94
message("Input mosaic copied to output raster: now, process boundary")
95
#browser()
96

    
97
# now, we SHOULD be able to modify outMosaic with new 'column values'.
98

    
99
# note: last north image and first south image rows are the same. 
100

    
101
     northBorderEdgeRow <- lastNorthSubImgRow
102
     southBorderEdgeRow <- firstSouthSubImgRow + 1
103

    
104
# The border 'edge' row is one row below the top of the south sub image
105

    
106
     southBrdrRowVals <- getValues(outMosaic,southBorderEdgeRow,1)
107
     northBrdrRowVals <- getValues(outMosaic,northBorderEdgeRow,1)
108
 
109
     brdrRowEdgeDelta <- southBrdrRowVals - northBrdrRowVals
110

    
111
# Process the mosaic 'column-by-column'
112
# For each column, extract the LAST 'numBlendRows cells adjacent to the bottomi
113
# (southernmost) edge of the ASTER component, next to the SRTM border. 
114
# Then, moving 'north' along the colum of values, add the current colunn's 'ledge' 
115
# offset, attenuated by the exponential decay factor for that row. 
116
# Thus, the SRTM/ASTER 'ledge' is attenuated, the impact is greatest at the boundary,
117
# least at the threshold, when the decay factor (empirically determined) goes to zero.
118
#
119
# Add Exponential Decay constant: empirically selected (try a few, ee if you like result
120
# Maybe plot it, see how many rows it takes for the 'effect' of SRTM goes to zero.
121
 
122
     dKc <- .045 # empirically determined constant, so that exp - 1 decay goes to zero in the desired number of image rows.
123
     numBlendRows <- 60
124
     begBlendRow <- (northBorderEdgeRow - numBlendRows) + 1
125
#     srows <- seq(southBorderEdgeRow,((southBorderEdgeRow + numBlendRows) - 1))
126
     srows <- seq((begBlendRow),northBorderEdgeRow)
127
# 
128
     for (curMosCol in 1 : ncol(mosaicLayers[[1]]))
129
     {
130
message(sprintf("transforming cur col: %d",curMosCol))
131
#browser()
132
#        deltaInc <- brdrRowEdgeDelta[curMosCol] / numBlendRows
133
        colVecCells <- cellFromRowCol(outMosaic,srows,curMosCol) 
134
        colVecValues <- getValuesBlock(outMosaic,row=begBlendRow,nrows=numBlendRows,col=curMosCol,ncols=1) 
135
        compareVec <- getValuesBlock(outMosaic,row=begBlendRow,nrows=numBlendRows+2,col=curMosCol,ncols=1) 
136
        curRowBoundaryOffset <- brdrRowEdgeDelta[curMosCol]
137

    
138
# note: we need to deal with NA values, which are possible. 
139

    
140
#        sumDeltaInc <- 0
141
#        for (iCtr in numBlendRows : 1) # this moves from 'south to north', decay increases with distance.
142
        for (iCtr in 1 : numBlendRows) # this moves from 'south to north', decay increases with distance.
143
#        for (iCtr in length(colVecValues) : 1) # this moves from 'south to north', decay increases with distance.
144
        {
145

    
146
# The idea: the closer to the border, the larger an increment we assign.
147
# in any case, increment the offset so that it is correct for any iupcoming non-NA column value
148

    
149
            if (!is.na(colVecValues[1]))
150
            {
151
#message(sprintf("cur col: %d Found a NON-NA vector",curMosCol))
152
#browser()
153
 
154
# exponential decay diminishes impact of the 'delta' edge adjustment with distance from the boundary.
155

    
156
               dK <- exp(-(iCtr * dKc))  
157
               colToMod <- (numBlendRows-iCtr)+1
158
#               qq <- as.integer(round(colVecValues[colToMod] + (curRowBoundaryOffset * dK)))
159
               colVecValues[colToMod] <- as.integer(round(colVecValues[colToMod] + (curRowBoundaryOffset * dK)))
160
#               colVecValues[iCtr] <- as.integer(round(colVecValues[iCtr] + (curRowBoundaryOffset * dK)))
161
            }
162
        }
163

    
164
# Insert this vector back into the mosaic: this technique adopted from raster PDF doc 'replacement' help. 
165
# Looks like: I have to specify the row index (vector) to insert the values into the mosaic. 
166
message("write columns to out mosaic")
167
#rowser()
168
        outMosaic[colVecCells] <- colVecValues
169

    
170
     }
171

    
172
# write the modified outMosaic values to the image file that we created.
173

    
174
message("Done with raster blending : hit key to update raster disk file..")
175
browser()
176
     writeRaster(outMosaic, sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
177
}
(4-4/18)