forked from mgimond/Spatial_pre2021
-
Notifications
You must be signed in to change notification settings - Fork 0
/
A05-Raster-Operations.Rmd
588 lines (429 loc) · 25.5 KB
/
A05-Raster-Operations.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
# Raster operations in R{-}
## Sample files for this exercise{-}
We'll first load spatial objects used in this exercise from a remote website: an elevation `raster` object, a bathymetry `raster` object and a continents `SpatialPolygonsDataFrame` vector layer.
```{r}
load(url("http://github.com/mgimond/Spatial/raw/master/Data/raster.RData"))
```
Both rasters cover the entire globe. Elevation below mean sea level are encoded as `0` in the elevation raster. Likewise, bathymetry values above mean sea level are encoded as `0`.
Note that most of the map algebra operations and functions covered in this tutorial are implemented using the `raster` package. See chapter 10 for a theoretical discussion of map algebra operations.
## Local operations and functions {-}
### Unary operations and functions (applied to single rasters) {-}
Most algebraic operations can be applied to rasters as they would with any vector element. For example, to convert all bathymetric values in `bath` (currently recorded as positive values) to *negative* values simply multiply the raster by `-1`.
```{r message=FALSE}
library(raster)
bath2 <- bath * (-1)
```
Another unary operation that can be applied to a raster is reclassification. In the following example, we will assign all `bath2` values that are less than zero a `1` and all zero values will remain unchanged. A simple way to do this is to apply a conditional statement.
```{r message=FALSE}
bath3 <- bath2 < 0
```
Let's look at the output. Note that all `0` pixels are coded as `FALSE` and all `1` pixels are coded as `TRUE`.
```{r, fig.height=3}
library(tmap)
tm_shape(bath3) + tm_raster(palette = "Greys") +
tm_legend(outside = TRUE, text.size = .8)
```
If a more elaborate form of reclassification is desired, you can use the `reclassify` function. In the following example, the raster object `bath` is reclassified to 4 unique values: 100, 500, 1000 and 11000 as follows:
Original depth values | Reclassified values
--------------------- | -------------------
0 - 100 | 100
101 - 500 | 500
501 - 1000 | 1000
1001 - 11000 | 11000
The first step is to create a plain matrix where the first and second columns list the starting and ending values of the range of input values that are to be reclassified, and where the third column lists the new raster cell values.
```{r}
m <- c(0, 100, 100, 100, 500, 500, 500,
1000, 1000, 1000, 11000, 11000)
m <- matrix(m, ncol=3, byrow = T)
m
```
```{r}
bath3 <- reclassify(bath, m, right = T)
```
The `right=T` parameter indicates that the intervals should be closed to the right (i.e. the second column of the reclassification matrix is inclusive).
```{r, fig.height=3}
tm_shape(bath3) + tm_raster(style="cat") + tm_legend(outside = TRUE, text.size = .8)
```
You can also assign `NA` (missing) values to pixels. For example, to assign `NA` values to cells that are equal to 100, type
```{r message=FALSE}
bath3[bath3 == 100] <- NA
```
The following chunk of code highlights all `NA` pixels in grey and labels them as missing.
```{r, fig.height=3}
tm_shape(bath3) + tm_raster(showNA=TRUE, colorNA="grey") +
tm_legend(outside = TRUE, text.size = .8)
```
### Binary operations and functions (where two rasters are used) {-}
In the following example, `elev` (elevation raster) is added to `bath` (bathymetry raster) to create a single elevation raster for the globe. Note that the bathymetric raster will need to be multiplied by `-1` to differentiate above mean sea level elevation from below mean sea level depth.
```{r message=FALSE, warning=FALSE, tidy=F}
elevation <- elev - bath
```
```{r, fig.height=3}
tm_shape(elevation) + tm_raster(palette="-RdBu",n=6) +
tm_legend(outside = TRUE, text.size = .8)
```
## Focal operations and functions {-}
Operations or functions applied focally to rasters involve user defined neighboring cells. For example, a cell output value can be the average of all 121 cells--an 11x11 kernel--centered on the cell whose value is being estimated (this is an example of a _smoothing function_).
```{r}
f1 <- focal(elevation, w=matrix(1,nrow=11,ncol=11) , fun=mean)
```
```{r, fig.height=3}
tm_shape(f1) + tm_raster(palette="-RdBu",n=6) +
tm_legend(outside = TRUE, text.size = .8)
```
By default edge cells are assigned a value of `NA`. This is because cells outside of the input raster extent have no value but may fall inside the kernel window when edge cells are computed. You can remedy this by passing the parameters `na.rm=TRUE` and `pad=TRUE` to the function.
```{r}
f1 <- focal(elevation, w=matrix(1,nrow=11,ncol=11) ,
fun=mean, pad=TRUE, na.rm = TRUE)
```
```{r, fig.height=3}
tm_shape(f1) + tm_raster(palette="-RdBu",n=6) +
tm_legend(outside = TRUE, text.size = .8)
```
But, you must be careful when making use of the `na.rm=TRUE`/`pad=TRUE` combination. In the above example, we are passing the values from the kernel window as *a vector element* to the `mean` function. So the number of valid (non NA) values being passed to the `mean` function does not matter since the function will sort out the weights based on the number of input values.
If, however, the smoothing function makes use of explicitly defined weights, then setting `pad=TRUE` will be problematic since the edge pixels will succumb to unbalanced weights. In the following example, we'll compare the output of three 3x3 smoothing operations as follows:
```{r}
# Using the mean function
f_mean <- focal(elevation, w=matrix(1,nrow=3,ncol=3), fun=mean, na.rm=TRUE, pad=TRUE)
# Using explicitly defined weights
f_wt_nopad <- focal(elevation, w=matrix(1/9,nrow=3,ncol=3), na.rm=TRUE, pad=FALSE)
f_wt_pad <- focal(elevation, w=matrix(1/9,nrow=3,ncol=3), na.rm=TRUE, pad=TRUE)
```
The following images are zoomed in on the upper left-hand set of pixels of each raster output. The `elevation` raster is the input raster.
```{r echo=FALSE, fig.width = 8, fig.height=1.5}
zoom1 <- extent(-180, -179, 89,90)
OP <- par(pty="s", mfrow=c(1,4), mar=c(0,0,3,0), bty="n")
plot(crop(elevation, zoom1), colNA="bisque", asp=1, axes=FALSE, legend=FALSE,
zlim=c(-4150, -2700), col=grey(1:9/12))
title("Elevation", line=1)
text(crop(elevation, zoom1), col="yellow")
grid(nx=3, ny=3, lty=1)
plot(crop(f_mean, zoom1), colNA="bisque", asp=1, axes=FALSE, legend=FALSE,
zlim=c(-4150, -2700), col=grey(1:9/12))
title("f_mean", line=1)
text(crop(f_mean, zoom1), col="yellow")
grid(nx=3, ny=3, lty=1)
plot(crop(f_wt_nopad, zoom1), colNA="bisque", asp=1, axes=FALSE, legend=FALSE,
zlim=c(-4150, -2700), col=grey(1:9/12))
title("f_wt_nopad", line=1)
text(crop(f_wt_nopad, zoom1), col="yellow")
grid(nx=3, ny=3, lty=1)
plot(crop(f_wt_pad, zoom1), colNA="bisque", asp=1, axes=FALSE, legend=FALSE,
zlim=c(-4150, -2700), col=grey(1:9/12))
title("f_wt_nopad", line=1)
text(crop(f_wt_pad, zoom1), col="yellow")
grid(nx=3, ny=3, lty=1)
par(OP)
```
The left-most image shows the original upper left-hand corner elevation pixel values. The second image shows the output using the mean function with the `pad=TRUE` option. The third image shows the output using the weighted kernel window with the `pad=FALSE` option (note the upper row of `NA` pixels highlighted in *bisque* color). The last image shows the output from the same weighted kernel window using the `pad=TRUE` option. This last outcome is not desirable. For example, the middle top pixel is computed from `1/9(-4113 -4113 -4112 -4107 -4104 -4103)`, which results in dividing the sum of *six* values by *nine*--hence the unbalanced weight effect. Note that we do not have that problem using the mean function.
You might have noticed the lack of edge effect issues along the western edge of the raster outputs. This is because the `focal` function will wrap the eastern edge of the raster to the western edge of that same raster if the input raster layer spans the entire globe (i.e from -180 ° to +180 °).
The _neighbors_ matrix (or _kernel_) that defines the moving window can be customized. For example if we wanted to compute the average of all 8 neighboring cells _excluding_ the central cell we could define the matrix as follows:
```{r eval=FALSE}
m <- matrix(c(1,1,1,1,0,1,1,1,1)/8,nrow = 3)
f2 <- focal(elevation, w=m, fun=sum)
```
More complicated kernels can be defined. In the following example, a [Sobel][2] filter (used for edge detection in image processing) is defined then applied to the raster layer `elevation`.
```{r fig.height=3}
Sobel <- matrix(c(-1,0,1,-2,0,2,-1,0,1) / 4, nrow=3)
f3 <- focal(elevation, w=Sobel, fun=sum)
tm_shape(f3) + tm_raster(palette="Greys") +
tm_legend(legend.show = FALSE)
```
## Zonal operations and functions {-}
A common zonal operation is the aggregation of cells. In the following example, raster layer `elevation` is aggregated to a 5x5 raster layer.
```{r fig.height=3}
z1 <- aggregate(elevation, fact=2, fun=mean, expand=TRUE)
tm_shape(z1) + tm_raster(palette="-RdBu",n=6) +
tm_legend(outside = TRUE, text.size = .8)
```
The image may not look much different from the original, but a look at the image properties will show a difference in pixel sizes.
```{r}
res(elevation)
res(z1)
```
`z1`'s pixel dimensions are half of `elevation`'s dimensions. You can reverse the process by using the `disaggregate` function which will split a cell into the desired number of subcells while assigning each one the same parent cell value.
Zonal operations can often involve two layers, one with the values to be aggregated, the other with the defined zones. In the next example, `elevation`'s cell values are averaged by zones defined by the `cont` polygon layer.
The following chunk computes the mean elevation value for each unique polygon in `cont`,
```{r}
cont.elev <- extract(elevation, cont, fun=mean, sp=TRUE)
```
The `sp=TRUE` parameter instructs the function to output a `SpatialPolygonsDataFrame` object (same as the input zonal object) instead of a standalone table (matrix). The output spatial object inherits the original attributes and adds a column called `layer` with the computed mean elevation values.
```{r}
cont.elev@data
```
We can now map the average elevation by continent.
```{r, fig.height=3}
tm_shape(cont.elev) + tm_polygons(col="layer") +
tm_legend(outside = TRUE, text.size = .8)
```
Many custom functions can be applied to `extract`. For example, to extract the maximum elevation value by continent, type:
```{r}
cont.elev <- extract(elevation, cont, fun=max, sp=TRUE)
```
As another example, we may wish to extract the number of pixels in each polygon using a customized function.
```{r}
cont.elev <- extract(elevation, cont, fun=function(x,...){length(x)}, sp=TRUE)
```
The `extract` function will also work with lines and point spatial objects. If you wish to compute the zonal statistics of a raster using another raster as zones instead of a vector layer, use the `zonal()` function instead.
## Global operations and functions {-}
Global operations and functions may make use of all input cells of a grid in the computation of an output cell value.
An example of a global function is the Euclidean distance function, `distance`, which computes the shortest distance between a pixel and a source (or destination) location. To demonstrate the `distance` function, we'll first create a new raster layer with two non-NA pixels.
```{r}
r1 <- raster(ncol=100, nrow=100, xmn=0, xmx=100, ymn=0, ymx=100)
r1[] <- NA # Assign NoData values to all pixels
r1[c(850, 5650)] <- 1 # Change the pixels #850 and #5650 to 1
```
```{r fig.height=3}
tm_shape(r1) + tm_raster(palette="red") +
tm_legend(outside = TRUE, text.size = .8)
```
Next, we'll compute a Euclidean distance raster from these two cells. The output extent will default to the input raster extent.
```{r}
r1.d <- distance(r1)
```
```{r fig.height=3}
tm_shape(r1.d) + tm_raster(palette = "Greens", style="order", title="Distance") +
tm_legend(outside = TRUE, text.size = .8) +
tm_shape(r1) + tm_raster(palette="red", title="Points")
```
You can also compute a distance raster using `SpatialPoints` objects or a simple x,y data table. In the following example, distances to points (25,30) and (87,80) are computed for each output cell. However, since we are working off of point objects (and not an existing raster as was the case in the previous example), we will need to create a blank raster layer which will define the extent of the Euclidean distance raster output.
```{r}
# Create a blank raster
r2 <- raster(ncol=100, nrow=100, xmn=0, xmx=100, ymn=0, ymx=100)
# Create a point layer
xy <- matrix(c(25,30,87,80),nrow=2, byrow=T)
p1 <- SpatialPoints(xy)
```
Now let's compute the Euclidean distance to these points using the `distanceFromPoints` function.
```{r}
r2.d <- distanceFromPoints(r2, p1)
```
Let's plot the resulting output.
```{r fig.height=3}
tm_shape(r2.d) + tm_raster(palette = "Greens", style="order") +
tm_legend(outside = TRUE, text.size = .8) +
tm_shape(p1) + tm_bubbles(col="red")
```
## Computing cumulative distances {-}
This exercise demonstrates how to use functions from the `gdistance` package to generate a cumulative distance raster. One objective will be to demonstrate the influence "adjacency cells" wields in the final results.
Load the `gdistance` package.
```{r}
library(gdistance)
```
First, we'll create a 100x100 raster and assign a value of `1` to each cell. The pixel value defines the cost (other than distance) in traversing that pixel. In this example, we'll assume that the cost is uniform across the entire extent.
```{r message=FALSE, warning=FALSE, tidy=F}
r <- raster(nrows=100,ncols=100,xmn=0,ymn=0,xmx=100,ymx=100)
r[] <- rep(1, ncell(r))
```
```{r echo=FALSE}
# Code block is ued internally to generate descriptive graphics.
library(ggplot2)
library(rasterVis)
library(reshape2)
# Plot raster subset (with cell values) and adjacency connections
plotsub = function(r, x1, x2, y1, y2){
r2 = crop(r, c(x1,x2,y1,y2))
gplot(r2) + geom_tile(aes(fill = value)) +
scale_colour_brewer(palette="Blues") +
coord_equal() +
geom_text(aes(label=sprintf("%2.1f", value)), colour='white',size=4)
}
plotsubconn = function(r, x1, x2, y1, y2, adj){
plotsub(r,x1,x2,y1,y2) +
geom_segment(data = adj, mapping=aes(x=x, y=y, xend=x+dx, yend=y+dy),
arrow = arrow(length=unit(0.1,"cm")),colour = "red") +
geom_point(data = adj, mapping=aes(x=x, y=y), size=12, shape=21, fill="yellow") +
geom_point(data = adj, mapping=aes(x=x+dx, y=y+dy), size=12, shape=21, fill="white") +
geom_text(aes(label=sprintf("%2.1f", value)), colour='black',size=4) +
theme(legend.position = "none")
}
# Plot raster subset (with cell values)
gp.text <- function(sg, nr=1, sfmt="%02.0f",x1,x2,y1,y2,src){
src = as.data.frame(t(src)); colnames(src)=c("x","y")
sg.names = names(sg)
sg = crop(sg, c(x1,x2,y1,y2))
names(sg) = sg.names # Renaming the layers is necessary if the brick object is cropped
r.df = as.data.frame(cbind(xyFromCell(sg,1:ncell(sg)),getValues(sg)))
plotData = melt(r.df, id=c("x","y"))
ggplot(aes(x = x, y = y), data = plotData, environment = environment()) +
geom_tile(aes(fill = value)) + facet_wrap(~ variable, nrow=nr) +
theme(strip.text.x = element_text(size=16, family="mono", face="bold"), legend.position = "none") +
geom_point(data = src, mapping=aes(x=x, y=y), size=18, shape=21, colour="yellow") +
coord_equal() + geom_text(aes(label=sprintf(fmt=sfmt, value)), colour='white',size=2)
}
# Plot rasters with extent limit
gp.notext <- function(sg, nr=1, x1,x2,y1,y2,lim,src){
src = as.data.frame(t(src)); colnames(src)=c("x","y")
sg.names = names(sg)
sg = crop(sg, c(x1,x2,y1,y2))
names(sg) = sg.names # Renaming the layers is necessary if the brick object is cropped
r.df = as.data.frame(cbind(xyFromCell(sg,1:ncell(sg)),getValues(sg)))
plotData = melt(r.df, id=c("x","y"))
ggplot(aes(x = x, y = y), data = plotData, environment = environment()) +
geom_tile(aes(fill = value)) + facet_wrap(~ variable, nrow=nr) +
theme(strip.text.x = element_text(size=16, family="mono", face="bold"), legend.position = "none") +
coord_equal() +
geom_tile(data = subset(plotData, value <= lim),aes(fill = value),fill="red") +
geom_point(data = src, mapping=aes(x=x, y=y), size=10, shape=21, fill="yellow")
}
```
```{r echo=FALSE}
r <- raster(nrows=100,ncols=100,xmn=0,ymn=0,xmx=100,ymx=100)
r[] <- rep(1, ncell(r))
# Define adjacency parameters (used only for visualizing neighbors)
rook <- data.frame(x=rep((nrow(r) + 1)/2, 4), y=rep((nrow(r) - 1)/2, 4),
dx=c(0,-1,1,0), dy=c(1,0,0,-1))
queen <- data.frame(x=rep((nrow(r) + 1)/2, 16), y=rep((nrow(r) - 1)/2, 16),
dx=c(-1,0,1,-1,1,-1,0,1), dy=c(1,1,1,0,0,-1,-1,-1))
queenknight <- data.frame(x=rep((nrow(r) + 1)/2, 8), y=rep((nrow(r) - 1)/2, 8),
dx=c(-1,1,-2,-1,0,1,2,-1,1,-2,-1,0,1,2,-1,1),
dy=c(2,2,1,1,1,1,1,0,0,-1,-1,-1,-1,-1,-2,-2))
bishop <- data.frame(x=rep((nrow(r) + 1)/2, 4), y=rep((nrow(r) - 1)/2, 4),
dx=c(-1,1,-1,1),
dy=c(1,1,-1,-1))
```
If you were to include traveling costs other than distance (such as elevation) you would assign those values to each cell instead of the constant value of `1`.
A **translation matrix** allows one to define a 'traversing' cost going from one cell to an adjacent cell. Since we are assuming there are no 'costs' (other than distance) in traversing from one cell to any adjacent cell we'll assign a value of 1, `function(x){1}`, to the translation between a cell and its adjacent cells (i.e. translation cost is uniform in all directions).
There are four different ways in which 'adjacency' can be defined using the `transition` function. These are showcased in the next four blocks of code.
In this example, adjacency is defined as a **four node** (vertical and horizontal) connection (i.e. a "rook" move).
```{r echo=FALSE, fig.height=3, fig.width=3, fig.align='left'}
plotsubconn(r,48,53,47,52,rook)
```
```{r message=FALSE, warning=FALSE, tidy=F}
h4 <- transition(r, transitionFunction = function(x){1}, directions = 4)
```
In this example, adjacency is defined as an **eight node** connection (i.e. a single cell "queen" move).
```{r echo=FALSE, fig.height=3, fig.width=3, fig.align='left'}
plotsubconn(r,48,53,47,52,queen)
```
```{r message=FALSE, warning=FALSE, tidy=F}
h8 <- transition(r, transitionFunction = function(x){1}, directions = 8)
```
In this example, adjacency is defined as a **sixteen node** connection (i.e. a single cell "queen" move combined with a "knight" move).
```{r echo=FALSE, fig.height=3, fig.width=3, fig.align='left'}
plotsubconn(r,48,53,47,52,queenknight)
```
```{r message=FALSE, warning=FALSE, tidy=F}
h16 <- transition(r, transitionFunction=function(x){1},16,symm=FALSE)
```
In this example, adjacency is defined as a **four node diagonal** connection (i.e. a single cell "bishop" move).
```{r echo=FALSE, fig.height=3, fig.width=3, fig.align='left'}
plotsubconn(r,48,53,47,52,bishop)
```
```{r message=FALSE, warning=FALSE, tidy=F}
hb <- transition(r, transitionFunction=function(x){1},"bishop",symm=FALSE)
```
The `transition` function treats all adjacent cells as being at an equal distance from the source cell across the entire raster. `geoCorrection` corrects for 'true' local distance. In essence, it's adding an additional cost to traversing from one cell to an adjacent cell (the original cost being defined using the `transition` function). The importance of applying this correction will be shown later.
Note: `geoCorrection` also corrects for distance distortions associated with data in a geographic coordinate system. To take advantage of this correction, make sure to define the raster layer's coordinate system using the `projection` function.
```{r echo=FALSE, message=FALSE, warning=FALSE}
h16.old <- h16 # Grab copy of h16 before applying correction for later use
```
```{r message=FALSE, warning=FALSE, tidy=F}
h4 <- geoCorrection(h4, scl=FALSE)
h8 <- geoCorrection(h8, scl=FALSE)
h16 <- geoCorrection(h16, scl=FALSE)
hb <- geoCorrection(hb, scl=FALSE)
```
In the "queen's" case, the diagonal neighbors are $\sqrt{2 x (CellWidth)^{2}}$ times the cell width distance from the source cell.
Next we will map the cumulative distance (`accCost`) from a central point (A) to all cells in the raster using the four different adjacency definitions.
```{r message=FALSE, warning=FALSE, tidy=F}
A <- c(50,50) # Location of source cell
h4.acc <- accCost(h4,A)
h8.acc <- accCost(h8,A)
h16.acc <- accCost(h16,A)
hb.acc <- accCost(hb,A)
```
If the `geoCorrection` function had not been applied in the previous steps, the cumulative distance between point location A and its neighboring adjacent cells would have been different. Note the difference in cumulative distance for the 16-direction case as shown in the next two figures.
Uncorrected (i.e. `geoCorrection` _not_ applied to h16):
```{r echo=FALSE, fig.height=3, fig.width=3, fig.align='left'}
h16.acc.nocorrect <- accCost(h16.old,A)
plotsubconn(h16.acc.nocorrect,48,53,47,52,queenknight)
```
Corrected (i.e. `geoCorrection` applied to h16):
```{r echo=FALSE, fig.height=3, fig.width=3, fig.align='left'}
plotsubconn(h16.acc,48,53,47,52,queenknight)
```
The "bishop" case offers a unique problem: only cells in the diagonal direction are identified as being adjacent. This leaves many undefined cells (labeled as `Inf`). We will change the `Inf` cells to `NA` cells.
```{r message=FALSE, warning=FALSE, tidy=F}
hb.acc[hb.acc == Inf] <- NA
```
Now let's compare a 7x7 subset (centered on point A) between the four different cumulative distance rasters.
```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height=6, fig.width=6, fig.align='left'}
sg <- brick(h4.acc,h8.acc,h16.acc,hb.acc)
names(sg) <- c("h4.acc","h8.acc","h16.acc","hb.acc")
gp.text(sg, nr=2, sfmt="%2.1f", x1=47,x2=54,y1=46,y2=53,src=c(50.5,49.5))
```
To highlight the differences between all four rasters, we will assign a red color to all cells that are within 20 cell units of point A.
```{r echo=FALSE, fig.height=6, fig.width=6, fig.align='left'}
sg <- brick(h4.acc,h8.acc,h16.acc,hb.acc)
names(sg) <- c("h4.acc","h8.acc","h16.acc","hb.acc")
gp.notext(sg, nr=2, x1=0,x2=100,y1=0,y2=100,lim=20,src=c(50.5,49.5))
```
```{r echo=FALSE, message=FALSE, warning=FALSE}
rb <- table(hb.acc[] <= 30)
r4 <- table(h4.acc[] <= 30)
r8 <- table(h8.acc[] <= 30)
r16 <- table(h16.acc[] <= 30)
Xr <- cbind (rb[2], r4[2], r8[2], r16[2])
```
It's obvious that the accuracy of the cumulative distance raster can be greatly influenced by how we define adjacent nodes. The number of red cells (i.e. area identified as being within a 20 units cumulative distance) ranges from `r min(Xr)` to `r max(Xr)` cells.
```{r echo=FALSE, fig.height=3}
barplot(Xr, names.arg=c("Bishop", "4 Node", "8 Node", "16 Node"),
col="bisque2",ylab="# of cells within a dist of 20 units")
```
### Working example {-}
In the following example, we will generate a raster layer with barriers (defined as `NA` cell values). The goal will be to identify all cells that fall within a 290 km traveling distance from the upper left-hand corner of the raster layer. Results between an 8-node and 16-node adjacency definition will be compared.
```{r fig.height=4, fig.align='left'}
# create an empty raster
r <- raster(nrows=300,ncols=150,xmn=0,ymn=0,xmx=150000, ymx=300000)
# Define a UTM projection (this sets map units to meters)
projection(r) = "+proj=utm +zone=19 +datum=NAD83"
# Each cell is assigned a value of 1
r[] <- rep(1, ncell(r))
# Generate 'baffles' by assigning NA to cells. Cells are identified by
# their index and not their coordinates.
# Baffles need to be 2 cells thick to prevent the 16-node
# case from "jumping" a one pixel thick NA cell.
a <- c(seq(3001,3100,1),seq(3151,3250,1))
a <- c(a, a+6000, a+12000, a+18000, a+24000, a+30000, a+36000)
a <- c(a , a+3050)
r[a] <- NA
# Let's check that the baffles are properly placed
tm_shape(r) + tm_raster(colorNA="red") +
tm_legend(legend.show=FALSE)
# Next, generate a transition matrix for the 8-node case and the 16-node case
h8 <- transition(r, transitionFunction = function(x){1}, directions = 8)
h16 <- transition(r, transitionFunction = function(x){1}, directions = 16)
# Now assign distance cost to the matrices.
h8 <- geoCorrection(h8)
h16 <- geoCorrection(h16)
# Define a point source.
A <- SpatialPoints(cbind(50,290000))
# Compute the cumulative cost raster
h8.acc <- accCost(h8, A)
h16.acc <- accCost(h16,A)
# Replace Inf with NA
h8.acc[h8.acc == Inf] <- NA
h16.acc[h16.acc == Inf] <- NA
```
Let's plot the results. Yellow cells will identify cumulative distances within 290 km.
```{r fig.height=4, fig.align='left'}
tm_shape(h8.acc) + tm_raster(n=2, style="fixed", breaks=c(0,290000,Inf)) +
tm_facets() + tm_shape(A) + tm_bubbles(col="red") +
tm_legend(outside = TRUE, text.size = .8)
tm_shape(h16.acc) + tm_raster(n=2, style="fixed", breaks=c(0,290000,Inf)) +
tm_facets() + tm_shape(A) + tm_bubbles(col="red") +
tm_legend(outside = TRUE, text.size = .8)
```
We can compute the difference between the 8-node and 16-node cumulative distance rasters:
```{r message=FALSE, warning=FALSE, tidy=F}
table(h8.acc[] <= 290000)
table(h16.acc[] <= 290000)
```
```{r echo=FALSE}
r8 <- table(h8.acc[] <= 290000)
r16 <- table(h16.acc[] <= 290000)
```
The number of cells identified as being within a 290 km cumulative distance of point A for the 8-node case is `r sprintf("%2.0f",r8[2])` whereas it's `r sprintf("%2.0f",r16[2])` for the 16-node case, a difference of `r sprintf("%2.1f%%",(1 - r8[2]/r16[2]) * 100)`.
[1]: http://en.wikipedia.org/wiki/Dana_Tomlin
[2]: http://en.wikipedia.org/wiki/Sobel_operator