Project

General

Profile

Download (7.1 KB) Statistics
| Branch: | Revision:
1
##################################################################################### 
2
# DEMBoundarySolutionsBlock - Block version, uses a BLOCK write to create output mosaic
3
# 
4
# Collection of functions that attempt to mitigate ASTER / SRTM the boundary edge problem
5
# 
6
# 
7
#
8
# Inputs: 
9
#   1) Blended Aster/SRTM DEM mosaic image.
10
#   2) 'Aster-Only' subimage corresponding to the NORTHernmost part (adjacent to the 60 degree Lat 
11
#     boundary) of the SRTM component of the mosaic.
12
#   3) 'Aster-Only' subimage corresponding to the SOUTHernmost part (adjacent to the 60 degree Lat 
13
#     boundary) of the ASTER component of the mosaic.
14
#   4) Other parameters, as they are needed
15
#
16
# To run:
17
#
18
#  1) start R
19
#  2) > source("DEMBoundarySolutions.r")
20
#  3) > BlendAsterDEMPixelsBlock()
21
#
22
# Author: Rick Reeves, NCEAS
23
# May 31, 2011
24

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

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

    
50
#     mosaicLayers <- stack("./AsterSrtmBoth3ArcSecSub.tif",
51
#                           "./AsterBdyTest3ArcSecFullGridAbove60.tif",
52
#                           "./AsterBdyTest3ArcSecFullGridBelow60.tif") # fix by exporting from ArcGIS
53

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

    
67
     northAsterEx <- extent(inAsterNorth)
68
     southAsterEx <- extent(inAsterSouth)
69
     
70
# Get the values from the mosaic for this extent
71

    
72
     southCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],southAsterEx)
73
     northCellsOfInterest <- cellsFromExtent(mosaicLayers[[1]],northAsterEx)
74

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

    
79
     firstNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[1])   
80
     lastNorthSubImgRow <- rowFromCell(mosaicLayers[[1]],northCellsOfInterest[(length(northCellsOfInterest) - 1)])
81
     northNrowsToProcess <- nrow(inAsterNorth)
82
  
83
     firstSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[1])   
84
     lastSouthSubImgRow <- rowFromCell(mosaicLayers[[1]],southCellsOfInterest[(length(southCellsOfInterest) - 1)])
85
     southNrowsToProcess <- nrow(inAsterSouth)
86

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

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

    
105
# note: last north image and first south image rows are the same. 
106

    
107
     northBorderEdgeRow <- lastNorthSubImgRow
108
     southBorderEdgeRow <- firstSouthSubImgRow + 1
109

    
110
# The border 'edge' row is one row below the top of the south sub image
111

    
112
     southBrdrRowVals <- getValues(outMosaic,southBorderEdgeRow,1)
113
     northBrdrRowVals <- getValues(outMosaic,northBorderEdgeRow,1)
114
 
115
     brdrRowEdgeDelta <- southBrdrRowVals - northBrdrRowVals
116

    
117
# Process the mosaic 'column-by-column'
118
# For first effort, apply straightforward linear ramp function to the northernmost SRTM image rows.
119
# In next iteration, we will modify the LAST 'numBlendRows' rows NORTH of the boundary with an 'extinction' function.
120
#
121
     numBlendRows <- 100
122
     srows <- seq(southBorderEdgeRow,((southBorderEdgeRow + numBlendRows) - 1))
123
# 
124
     for (curMosCol in 1 : ncol(mosaicLayers[[1]]))
125
     {
126
message(sprintf("transforming cur col: %d",curMosCol))
127
#browser()
128
        deltaInc <- brdrRowEdgeDelta[curMosCol] / numBlendRows
129
        colVecCells <- cellFromRowCol(outMosaic,srows,curMosCol) 
130
        colVecValues <- getValuesBlock(outMosaic,row=southBorderEdgeRow,nrows=numBlendRows,col=curMosCol,ncols=1) 
131

    
132
# note: we need to deal with NA values, which are possible. 
133

    
134
        sumDeltaInc <- 0
135
        for (iCtr in length(colVecValues) : 1)
136
        {
137

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

    
141
            sumDeltaInc <- sumDeltaInc + deltaInc
142

    
143
            if (!is.na(colVecValues[1]))
144
            {
145
#message(sprintf("cur col: %d Found a NON-NA vector",curMosCol))
146
#browser()
147
               colVecValues[iCtr] <- as.integer(round(colVecValues[iCtr] - sumDeltaInc)) 
148
            }
149
        }
150

    
151
# Insert this vector back into the mosaic: this technique adopted from raster PDF doc 'replacement' help. 
152
# Looks like: I have to specify the row index (vector) to insert the values into the mosaic. 
153
#message("write columns to out mosaic")
154
#browser()
155
        outMosaic[colVecCells] <- colVecValues
156
#        outMosaic[srows,curMosCol] <- colVecValues #  this did not seem to work - june 1
157

    
158
     }
159
#message(sprintf("cur col: %d about to write to outMosaic",curMosCol))
160
#browser()
161

    
162
# write the modified outMosaic values to the image file that we created.
163

    
164
message("Done with raster blending : hit key to update raster disk file..")
165
browser()
166
     writeRaster(outMosaic, sFixedRasterFile, datatype="INT2S", format="GTiff",overwrite=TRUE)
167
}
(5-5/19)