class: right, bottom

## Raster data in R

##### Nikhil Kaza
##### Department of City & Regional Planning <br /> University of North Carolina at Chapel Hill
###### updated: 2018-08-03

---
# Fields vs. Objects

This is an age old debate as to which is a better representation of reality.

Objects: Discrete, sharply defined boundaries, have distinct attributes (e.g. buildings, roads, parcels, census tracts)
  - Need to keep track of topology and spatial relations among them

Fields: Something that varies continuously over space. The discretisation is an artefact of the data storage and representation (e.g. temperature, ground level, urbanicity)
  - Topological relationships are embedded 
 
Usually fields are represented by rasters. Rasters could be thought of as a 2D-matrix of values, with row and column indices representing the location information.
---

# Raster data in urban analytics

.pull-left[
- Satellite images

- Ground level pictures

- Airborne sensors

- Radar
]

.pull-right[
![](./figs/examplerasters.png)
[Image Credit: ESRI](https://developers.arcgis.com/net/latest/uwp/guide/add-raster-data.htm)
]

---

# Raster package

In this course we will use Robert Hijmans' excellent Raster package. This package handles raster data that does not fit into memory.

```r
if(!require(raster)){
    install.packages("raster")
*    library(raster) 
}

(rast <- raster())
```

```
## class       : RasterLayer 
## dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
```

---

```r
res(rast) <- 20
rast
```

```
## class       : RasterLayer 
## dimensions  : 9, 18, 162  (nrow, ncol, ncell)
## resolution  : 20, 20  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
```

```r
ncol(rast) <- 30
 rast
```

```
## class       : RasterLayer 
## dimensions  : 9, 30, 270  (nrow, ncol, ncell)
## resolution  : 12, 20  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
```
---

# External Files

```r
rast <- raster(system.file("external/test.grd", package="raster"))
library(rasterVis)
levelplot(rast, xlab=NULL, ylab=NULL, scales=list(draw=FALSE), margin=FALSE,  colorkey= list(labels=list(cex=2.5)))
```

![](figs/unnamed-chunk-3-1.png)
See [https://oscarperpinan.github.io/rastervis/#levelplot](https://oscarperpinan.github.io/rastervis/#levelplot)
---

# Rasters can be stacked

```r
rast2 <- rast * runif(ncell(rast)) 
(rastS <- stack(rast, rast2))
```

```
## class       : RasterStack 
## dimensions  : 115, 80, 9200, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 40, 40  (x, y)
## extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:28992 +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs 
## names       :       test.1,       test.2 
## min values  : 128.43400574,   0.04142424 
## max values  :      1805.78,      1583.89
```
---

# or bricked

```r
rastB <- brick(rast, rast2)
names(rastB) <- c('Original', 'Modified')
levelplot(rastB, xlab=NULL, ylab=NULL, scales=list(draw=FALSE), margin=FALSE,  colorkey= list(labels=list(cex=2.5)))
```

![](figs/unnamed-chunk-5-1.png)
---

# Cell by cell algebra

Standard mathematical operators +,- etc. or logical operators such as `$\ge$`, `$\max$` etc.. work cell by cell between multiple Raster* objects or between Raster and numbers

```r
rastB[[1]] <- rast * 10
(rastB)
```

```
## class       : RasterBrick 
## dimensions  : 115, 80, 9200, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 40, 40  (x, y)
## extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:28992 +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs 
## data source : in memory
## names       :         test,     Modified 
## min values  : 1.284340e+03, 4.142424e-02 
## max values  :     18057.80,      1583.89
```

---
# Cell by cell algebra
More complicated algebra is possible

```r
(rastB[[1]]/100 >= rastB[[2]]) %>%
  levelplot(xlab=NULL, ylab=NULL, scales=list(draw=FALSE), margin=FALSE,  colorkey= list(labels=list(cex=2.5)))
```

---
# Moving window

```r
# 3x3 mean filter
(rast3 <- focal(rast, w=matrix(1/9,nrow=3,ncol=3)) )
```

```
## class       : RasterLayer 
## dimensions  : 115, 80, 9200  (nrow, ncol, ncell)
## resolution  : 40, 40  (x, y)
## extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:28992 +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : 168.611, 1372.349  (min, max)
```
---

# Moving window

```r
# Local maximum in 5X5 neighborhood
rast3 <- focal(rast, w=matrix(1, nrow=5,ncol=5), fun=max, na.rm=TRUE)
res(rast3)
```

```
## [1] 40 40
```

```r
levelplot(rast3,xlab=NULL, ylab=NULL, scales=list(draw=FALSE), margin=FALSE,  colorkey= list(labels=list(cex=2.5)))
```

<img src="figs/unnamed-chunk-9-1.png" width="70%" />
---

# Reduce resolution

```r
rast3 <- rast %>% aggregate(fact=5,fun=mean) 
res(rast3)
```

```
## [1] 200 200
```

```r
levelplot(rast3, xlab=NULL, ylab=NULL, scales=list(draw=FALSE), margin=FALSE,  colorkey= list(labels=list(cex=2.5)))
```

---

# Categorical raster

```r
rastCat <- rast %>% cut(breaks=c(0,300,800,1100, 1900))
rastCat <- ratify(rastCat)
rat <- levels(rastCat)[[1]]
rat$class <- c('Low', 'Mid', 'High', 'V. High')
```
<img src="figs/unnamed-chunk-12-1.png" width="60%" />

---

# Vector/Raster Operations

```r
cds1 <- structure(c(179169.418312112, 179776.022915626, 179344.40040928, 
179064.429053812, 178737.795805766, 179076.094526956, 178726.130332621, 
179169.418312112, 330830.748571234, 330527.446269478, 330049.161870553, 
329955.838085397, 330037.496397409, 330247.47491401, 330387.460591744, 
330830.748571234), .Dim = c(8L, 2L))

cds2 <- structure(c(180440.954884862, 180522.613196873, 179857.681227637, 
179951.005012793, 180440.954884862, 332417.252918885, 331717.324530216, 
331554.007906193, 331962.29946625, 332417.252918885), .Dim = c(5L, 
2L))

polys <- SpatialPolygonsDataFrame(SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                                  Polygons(list(Polygon(cds2)), 2))),data.frame(ID=c(1,2)))

polys
```

```
## class       : SpatialPolygonsDataFrame 
## features    : 2 
## extent      : 178726.1, 180522.6, 329955.8, 332417.3  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 1
## names       : ID 
## min values  :  1 
## max values  :  2
```
---

# Vector/Raster Operations

```r
p + layer(sp.polygons(polys, lwd=4, col='blue'))
```

![](figs/unnamed-chunk-14-1.png)

---
# Count percentage of cells

```r
library(tabularaster)
cn <- cellnumbers(rastCat, polys)
library(dplyr)
(cn %>% mutate(v = raster::extract(rastCat, cell_)) %>% 
    group_by(object_, v) %>% 
    summarize(count = n()) %>% 
  mutate(v.pct = count / sum(count))
  )
```

```
## # A tibble: 4 x 4
## # Groups:   object_ [2]
##   object_     v count v.pct
##     <int> <dbl> <int> <dbl>
## 1       1     1   114 0.368
## 2       1     2   196 0.632
## 3       2     1    49 0.249
## 4       2     2   148 0.751
```
---

# Least cost path between two points

```r
library(gdistance)
# Construct cells that are `impasssable`
cost <- calc(rast, fun=function(x)  ifelse(x > 400 & x<900, NA,  x)) #toy example
rclmat <- c(0,200, 20, 200, 400, 30, 900,1900, 40) %>% matrix(ncol=3, byrow = TRUE)
cost <- reclassify(cost, rclmat  )
A <- c(179706.2, 330570.0)
B <- c(180100.9, 331074.3)
plot(cost)
```

---

```r
# Least cost path between two points

conductance <- transition(cost, function(x) 1/mean(x), 8) # Conductance matrix
AtoB <- shortestPath(conductance, A, B, output="SpatialLines")
plot(cost)
lines(AtoB, col="red", lwd=2)
```

---
class: center, middle

Slides created via the R package [**xaringan**](https://github.com/yihui/xaringan).

The chakra comes from [remark.js](https://remarkjs.com), [**knitr**](http://yihui.name/knitr), and [R Markdown](https://rmarkdown.rstudio.com).