Tasks

Task 1 Write some code to determine the mean value and standard deviation of median income variable MedInc from the georgia data. Hints: the mean and sd functions can help you!

Task 2 Write some code to subset and plot the counties in the 25 percentile of median income from the georgia data. Hint: the quantile function can help you and you may want to assign the result to an R object that you can use to determine a threshold with which to subset the data as was done with index above, using a logical statement.

Task 3

Finally, go back to Session 3 and examine how the tb$rural variable was created for the tb tibble object. The same operations can be used to on the data.frame of sp objects. You should create a variable called rural that separates the counties of georgia into Rural and Non-Rural classes and determine how many counties are Rural and Non-Rural using this definition. Hint the table function can help here.

A quick look

The code below loads a dataset that is included in the sf package. It is of North Carolina (see https://cran.r-project.org/web/packages/spdep/vignettes/sids.pdf) and is loaded using the st_read function. This is a polygon spatial dataset…

# read the data
nc <- st_read(system.file("shape/nc.shp", package="sf"))

## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.3/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID):    4267
## proj4string:    +proj=longlat +datum=NAD27 +no_defs

which can of course be mapped

# do a quick tmap of the data
qtm(nc)

Quick tmap maps of North Carolina.

# and others mapping attributes
qtm(nc, fill = c("BIR74", "NWBIR79"))

Quick tmap maps of North Carolina.

Note that here we are using the qtm (quick tmap) function from the tmap package and not the plot function as above but also note that tmap can work with both sp and sf data formats. We can examine the data in the same way as a-spatial data:

dim(nc)
head(nc)
class(nc)

And again, you can see that nc has 100, 15 rows and columns respectively, that the class of nc indicate that it is sf class as well as data.frame. For data.frame related formats, R prints out all records and all columns. Try entering:

nc

Again you can generate and examine summaries of the nc data:

summary(nc)

Finally, it is also possibly to stop the legend and the body of the map from overlapping - although this can leave less space for the map itself:

qtm(nc, fill = "BIR74") +
  tm_layout(legend.outside = TRUE) 

qtm(nc, fill = "NWBIR79") +
  tm_layout(legend.outside = TRUE) 

sf vignettes

The sf package has a number of vignettes or tutorials that you could explore. These include an overview of the format, reading and writing from and to sf formats including conversions to and from sp and sf, and some illustrations of how sf objects can be manipulated. The code below will create a new window with a list of sf vignettes:

library(sf)
vignette(package = "sf")

To display a specific vignette topic, this can be called using the vignette function. A full description of the sf format can be found in thesf1` vignette

vignette("sf1", package = "sf")

Vignettes are an important part of R packages. They provide additional explanations of the package functionality to that found in the example code at the end of a help page. They can be accessed using the vignette function or through the R help. The sf1 vignette could also be accessed via the package help index: enter help(sf), then navigate to the index through the link at the bottom of the overview page and then click on the User guides, package vignettes and other documentation link.

Areas

It is possible to generate an outline of the area using the st_union() function in sf. An alternative for sp' is the gUnaryUnion() function in the rgeos package loaded with GISTools. The manipulation of spatial data using overlay, union and intersection functions is covered in more depth in Chapter 5 of Comber and Brunsdon (2015).

# do a merge
g <- st_union(georgia_sf)
# for sp
#georgia.outline <- gUnaryUnion(georgia, id = NULL)

Now run the code snippets below and inspect the results. You should see how the tmap functions are added as a series of layers to the map in a similar way to ggplot.

# plot the spatial layers
tm_shape(georgia_sf) +
  tm_fill("indianred") 

Add the border:

tm_shape(georgia_sf) +
  tm_fill("indianred") +
  tm_borders(lty = "dashed", col = "khaki")

Adding some styling:

tm_shape(georgia_sf) +
  tm_fill("indianred") +
  tm_borders(lty = "dashed", col = "khaki") +
  tm_style_natural(bg.color = "grey75") 

Adding the outline:

tm_shape(georgia_sf) +
  tm_fill("indianred") +
  tm_borders(lty = "dashed", col = "khaki") +
  tm_style_natural(bg.color = "grey75") +
  # now add the outline
  tm_shape(g) +
  tm_borders(lwd = 3) 

And finally putting it all together:

tm_shape(georgia_sf) +
  tm_fill("indianred") +
  tm_borders(lty = "dashed", col = "khaki") +
  tm_style_natural(bg.color = "grey75") +
  # now add the outline
  tm_shape(g) +
  tm_borders(lwd = 3) +
  tm_layout(title = "The State of Georgia", 
    title.size = 1, title.position = c(0.55, "top"))

So what you can see in the above code are two sets of tmap plot commands: the first set plots the georgia_sf dataset, specifying a dashed blue line to show the county boundaries, a red fill colour for the county and a map background colour of wheat. The second set adds the outline created by the union operation with a thicker line width before the title and subtitle are added. It is also possible to plot multiple different maps from different datasets together but this requires a bit more control over the tmap parameters. Note the way that the aspect of the second plot is used. The value was determined through trial and error. The code below assigns the map to a variable, and then a second set of functions, is used to manipulate these in a plot window. You may need to install the grid package:

# 1st plot of georgia
t1 <- tm_shape(georgia_sf) +
  tm_fill("coral") +
  tm_borders() +
  tm_style_natural(bg.color = "grey55") 
# 2nd plot of georgia2
t2 <- tm_shape(georgia2) +
  tm_fill("orange") +
  tm_borders() +
  tm_style_natural(bg.color = "grey75")+
  # the asp paramter controls aspect
  # this is makes the 2nd plot align
  tm_layout(asp = 0.86)

Now you can specify the layout of the combined map plot:

# open a new plot page
grid.newpage()
# set up the layout
pushViewport(viewport(layout=grid.layout(1,2)))
# plot using the print command
print(t1, vp=viewport(layout.pos.col = 1, height = 5))
print(t2, vp=viewport(layout.pos.col = 2, height = 5))

Thus different plot parameters can be used for different subsets of the data such that they are plotted in ways that are different from the default. Sometimes we would like to label the features in our maps. Have a look at the names of the counties in the georgia_sf dataset. These are held in the 13th attribute column and names(georgia_sf) will return a list of the names of all attributes:

data.frame(georgia_sf)[,13]

It would be useful to display these on the map and this can be done using the tm_text function in the maptools package that is loaded with tmap. The result is shown in the figure below.

tm_shape(georgia_sf) +
  tm_fill("white") +
  tm_borders() +
  tm_text("Name", size = 0.3) +
  tm_layout(frame = FALSE)

And we can subset the data as with sp format. The code below subsets the counties of Jefferson, Jenkins, Johnson, Washington, Glascock, Emanuel, Candler, Bulloch, Screven, Richmond and Burke.

# the county indices below were extracted from the data.frame
index <- c(81, 82, 83, 150, 62, 53, 21, 16, 124, 121, 17)
georgia_sf.sub <- georgia_sf[index,]

The notation is the same as for `sp object earlier, re-enforcing the fact that it is possible to select individual areas or polygons from spatial datasets using the bracket notation as used in matrices and vectors.

The subset can be plotted:

tm_shape(georgia_sf.sub) +
  tm_fill("gold1") +
  tm_borders("grey") +
  tm_text("Name", size = 1) +
  # add the outline
  tm_shape(g) +
  tm_borders(lwd = 2) +
  # specify some layout parameters
  tm_layout(frame = FALSE, title = "A subset of Georgia", 
    title.size = 2, title.position = c(0., "bottom"))

Finally, we can bring the different spatial data that has been created together in a single map. In the code below you should note how the different tm_shape, tm_fill etc functions are used to set up each layer of the map and that tmap determines the map extent from the layers:

# the 1st layer
tm_shape(georgia_sf) +
  tm_fill("white") +
  tm_borders("grey", lwd = 0.5) +
  # the 2nd layer
  tm_shape(g) +
  tm_borders(lwd = 2) +
  # the 3rd layer
  tm_shape(georgia_sf.sub) +
  tm_fill("lightblue") +
  tm_borders() +
  # specify some layout parameters
  tm_layout(frame = T, title = "Georgia with a subset of counties", 
    title.size = 1, title.position = c(0.02, "bottom"))

Points

Here we will briefly deal with mapping points using tmap. The code below creates some points from the coordinates of the georgia_sf polygons centroids (ignore the warnings!)

pts_sf <- st_centroid(georgia_sf)

## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data

These can be plotted in the same way as before but rather than using tm_polygons, tm_fill or tm_borders the tm_bubbles or tm_dots aesthetics can be specified:

# the 1st layer
tm_shape(georgia_sf) +
  tm_fill("olivedrab4") +
  tm_borders("grey", lwd = 1) +
  # the point layer
  tm_shape(pts_sf) +
  tm_dots(size = 1, alpha = 0.5)

Note how in the code above the tm_shape function was invoked for each layer.

Areas

It is also possible to plot multiple attributes simply by passing them to the tm_fill function. The plot window may need to be expanded to accommodate the legend:

tm_shape(georgia_sf) +
  tm_fill(c("PctRural", "PctBlack")) +
  tm_borders() +
  tm_layout(legend.format = list(digits = 0), 
            legend.position = c("right", "top"),
            legend.text.size = 0.5,
            legend.title.size = 0.8)

Note the number of plot parameters that are set by tm_legend - it is worthwhile exploring the help for this:

?tm_layout

You could also try legend.position = c("centre","bottom")). It looks bad! Further documentation ontm_layout` can be found at

Try repeating the tmap code above using a different variable such as PctEld. What happens? You will see that tmap automatically determines the number of classes to be included and the class intervals or breaks. Finally a colour shading scheme is automatically allocated to the map and the legend is included in the map. All of these, and many other of the default mapping settings that tmmap uses can be controlled and modified.

The tmap package has inbuilt tools change the way the class interval boundaries are computed. As a default, they are based on equal-sized intervals of the attribute being mapped, and different palette styles are available. Have a look at the help for tm_polygons and you will see that a number of different plotting styles are available. You should explore these.

The class intervals can be changed to quantiles or any other range of intervals using the breaks parameter. For example, the code below produces 2 maps with equal intervals (left), intervals based on k-means and with quantiles as (right).

However, the maps are not printed out. Instead they are assigned to p1 andp2. These can be called directly - try entering p2 at the console or used to create a multiple map figure using the pushViewport function in the grid package:

# with equal interval
p1 <-tm_shape(georgia_sf) +
  tm_polygons("PctBlack", title = "% Black", palette = "Greys") +
  tm_layout(legend.title.size = 1)
# with a style = kmeans
p2 <- tm_shape(georgia_sf) +
  tm_polygons("PctRural", title = "% Rural", palette = "Oranges",
    style = "kmeans") +
  tm_layout(legend.title.size = 1)

A multiple map figure can be created using the pushViewport function as below

# Multiple plots using the grid package
grid.newpage()
# set up the layout (rows, columns)
pushViewport(viewport(layout=grid.layout(1,2)))
# plot using he print command
print(p1, vp=viewport(layout.pos.col = 1, height = 5))
print(p2, vp=viewport(layout.pos.col = 2, height = 5))

Finally the class intervals can be controlled by manually specifying the breaks parameter:

tm_shape(georgia_sf) +
  tm_polygons("PctPov", breaks=seq(0, 100, by=25))

and this can be done in many different ways:

tm_shape(georgia_sf) +
  tm_polygons("PctPov", breaks=c(1,5,10, 20,30,40))

It is also possible to alter the colours used in a shading scheme. The default colour scheme uses increasing intensities of yellow red. You will have noticed the use of the different palettes in the above. The RColorBrewer package has a number of graduated lists of colours and is automatically loaded with both tmap and GISTools. The palettes designed by Cynthia Brewer and they are intended to optimise the perceptual difference between each shade in the palette, so that visually each shading colour is distinct (sometimes even when converted to a greyscale). The palettes available in this package are displayed with the command:

display.brewer.all()

This displays the various colour palettes and their names in a plot window. To generate a list of colours from one of these palettes, for example enter the following:

brewer.pal(5,'Blues')

## [1] "#EFF3FF" "#BDD7E7" "#6BAED6" "#3182BD" "#08519C"

This is a list of colour codes used by R to specify the palette. The arguments to brewer.pal specify that a five-stage palette based on shades of blue is required. The output of brewer.pal can be fed into tmap to give alternative colours in shading schemes. For example, enter the code below and a choropleth map shaded in green is displayed with its legend. The palette argument in tm_polygons specifies the new colours in the shading scheme.

tm_shape(georgia_sf) +
  tm_polygons("MedInc", title = "Median Income", palette = "Reds") +
  tm_layout(legend.title.size = 1, frame = F)

Note that the same map would be produced if the tm_fill function was used instead of tm_polygons in this case, however without a tm_borders function, the census block outlines are not plotted:

tm_shape(georgia_sf) +
  tm_fill("PctEld", title = "Median Income", palette = "Blues") +
  tm_layout(legend.title.size = 1, frame = F, legend.format = list(digits = 1))

It is also possible to display a histogram of the distribution of the variable or attribute being mapped using the legend.hist parameter. This is a very useful for choropleth mapping as it gives a distribution of the attributes being examined. Bringing this all together allows you to create a map that with a number of refinements as in the figure below. Note for example the minus sign before the palette parameter is used to reverse the palette order and the various parameters passed to the tm_layout function.

tm_shape(georgia_sf) +
  tm_polygons("PctEld", title = "Percentage Elderly", palette = "-GnBu",
              style = "kmeans",
              legend.hist = T) +
  tm_scale_bar(width = 0.22, position = c("right", "top")) +
  tm_compass()+
  tm_layout(frame = F, title = "Georgia", 
            title.size = 2, 
            #title.position = c(0.55, "top"), 
            legend.hist.size = 0.5, 
            legend.outside = T) 

New variables and attributes can be computed and mapped attributes values on the fly in tmap. The code below plots population density using the convert2density function applied to the TotPop90 attribute.

tm_shape(georgia_sf) +
    tm_fill(col="TotPop90", convert2density=TRUE,
        style="kmeans", title = expression("Population (per " * km^2 * ")"), 
        legend.hist=F, id="name") +
    tm_borders("grey25", alpha=.5) + 
  tm_format_NLD_wide()

The convert2density function automatically converts the projection units (in this case degrees of latitude and longitude) to a projection in metres and then determines areal density in km². You can check this by creating your own population density values, and examining the explanations of how the functions operate in the help pages for the functions used such as st_area.

Final note: here we are just using variables that are already included with the R spatial data objects. Of course you can create new variables in your analysis and assign them to the an sp or sf object.

The code below assigns some random numbers to a variable and generates a choropleth map:

georgia_sf$random <- rnorm(nrow(georgia_sf))
tm_shape(georgia_sf) +
    tm_polygons(col="random", style="kmeans", 
            title = "Random", palette = "RdBu") +
    tm_format_NLD_wide()

In summary, the tm_fill and tm_polygons functions in the tmap package generate choropleth maps of attributes held in SpatialPolygonsDataFrame (sp) or Simple Feature (sf) data objects. They automatically shade the variables using equal intervals. The intervals and the palettes can both be adjusted. It is instructive to examine the plotting functions and the way they operate. You can see the code in any R function by simply calling it without any arguments or brackets:

tm_polygons

The function code detail is displayed in the R console window. You will see that it takes a number of arguments and a number of default parameters. In addition to using the R help system to understand functions, examining functions in this way can also provide you with insight into their operation.

Points

Again we will briefly deal with mapping values or attributes for points using tmap. There are 2 basic options: change the size or the shading of the points according to the attribute value. We can use the pts_sf created earlier to illustrate these approaches.

First by size using tm_bubbles:

# the 1st layer
tm_shape(georgia_sf) +
  tm_fill("olivedrab4") +
  tm_borders("grey", lwd = 1) +
  # the points layer
  tm_shape(pts_sf) +
  tm_bubbles("PctBlack", title.size = "% Black", col = "gold")+
  tm_format_NLD()

Second by shade using tm_dots:

# the 1st layer
tm_shape(georgia_sf) +
  tm_fill("grey70") +
  tm_borders("gold", lwd = 1) +
  # the points layer
  tm_shape(pts_sf) +
  tm_dots("PctFB", size = 0.7, title = "% Foreign", shape = 15) +
  tm_layout(legend.outside = T, frame = F)