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

Raster package

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

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

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

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

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)))

See https://oscarperpinan.github.io/rastervis/#levelplot

Rasters can be stacked

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

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)))

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

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

(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

# 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

# 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

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

Reduce resolution

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

## [1] 200 200

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

Categorical raster

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')

Vector/Raster Operations

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

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

Count percentage of cells

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

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)

# 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)