diff --git a/01-introduction.Rmd b/01-introduction.Rmd
index 374ea96d6..848e66407 100644
--- a/01-introduction.Rmd
+++ b/01-introduction.Rmd
@@ -73,10 +73,10 @@ class(world)
 ```
 
 Let's look the first 2 rows and 3 columns of this object.
-The output shows 2 major differences compared with a regular `data.frame`: the inclusion of additional geographical data (`geometry type`, `dimension`, `bbox` and CRS information - `espf (SRID)`, `proj4string`), and the presence of a 4^th^ `geometry` column:
+The output shows 2 major differences compared with a regular `data.frame`: the inclusion of additional geographical data (`geometry type`, `dimension`, `bbox` and CRS information - `espf (SRID)`, `proj4string`), and the presence of a 4^th^ `geometry` column (result not shown):
 
-```{r}
-summary(world[1:3])
+```{r, results='hide'}
+summary(world)
 ```
 
 All this may seem rather complex, especially for a class system that is supposed to be simple.
diff --git a/03-attr.Rmd b/03-attr.Rmd
index 50aa30a42..845cbb1cb 100644
--- a/03-attr.Rmd
+++ b/03-attr.Rmd
@@ -119,11 +119,11 @@ world = world %>%
         mutate(pop_density = population / area)
 ```
 
-```{r}
-# data summary
+```{r, results='hide'}
+# data summary (not shown)
 summary(world)
 
-# data summary by groups
+# data summary by groups (not shown)
 world_continents = world %>% 
         group_by(continent) %>% 
         summarise(continent_pop = sum(population), country_n = n())
diff --git a/04-basic-map.Rmd b/04-basic-map.Rmd
index 66befb387..af8bf4a4f 100644
--- a/04-basic-map.Rmd
+++ b/04-basic-map.Rmd
@@ -12,10 +12,80 @@ plot(world)
 plot(world["population"])
 ```
 
+As with **sp**, you can add layers to your maps created with `plot()`, with the argument `add = TRUE`.^[In
+fact, when you `plot()` an **sf** object, R is calling `sf:::plot.sf()` behind the scenes.
+`plot()` is a generic method that behaves differently depending on the class of object being plotted.]
+However, this only works if the initial plot has only 1 layer (result not shown):
+
+```{r, fig.show='hide'}
+plot(world["pop_density"])
+africa = world_continents[1,]
+plot(africa, add = TRUE, col = "red")
+```
+
+This can be very useful when quickly checking the geographic correspondence between two or more layers.
+These plots work well for gaining a quick understanding of the data with few lines of code.
+For more advanced map making we recommend using a dedicated visualisation package such as **tmap**.
 
 <!-- 
 - plot() function 
 - map export 
 -->
 
+## Challenge
+
+Using **sf**'s `plot()` command, create a map of Nigeria in context, like the one presented in figure \@ref(fig:nigeria). 
+
+- Hint: this used the `lwd`, `main` and `col` arguments of `plot()`. 
+- Bonus: make the country boundaries a dotted grey line.
+- Hint: `border` is an additional argument of `plot()` for **sf** objects.
+
+```{r nigeria, warning=FALSE, echo=FALSE, fig.cap="Map of Nigeria in context illustrating sf's plotting capabilities"}
+nigeria = filter(world, name == "Nigeria")
+bb_africa = st_bbox(africa)
+plot(africa[2], col = "white", lwd = 3, main = "Nigeria in context", border = "lightgrey")
+# plot(world, lty = 3, add = TRUE, border = "grey")
+plot(world, add = TRUE, border = "grey")
+plot(nigeria, col = "yellow", add = TRUE, border = "darkgrey")
+ncentre = st_centroid(nigeria)
+ncentre_num = st_coordinates(ncentre)
+text(x = ncentre_num[1], y = ncentre_num[2], labels = "Nigeria")
+```
+
+# Further work
+
+**sf** makes R data objects more closely alligned to the data model used in GDAL and GEOS, in theory making spatial data operations faster.
+The work here provides a taster of the way that **sf** operates but there is much more to learn.
+There is a wealth of information that is available in the package's vignettes: these are highly recommended.
+
+As a final exercise, we'll see how to do a spatial overlay in **sf** by first converting the countries of the world into centroids and then subsetting those in Africa:
+
+```{r, out.width="50%", fig.cap="Centroids in Africa"}
+world_centroids = st_centroid(world)
+plot(world_centroids[1])
+africa_centroids = world_centroids[africa,]
+plot(africa_centroids, add = TRUE, cex = 2)
+```
+
+Note: another way of acheiving the same result is with a GEOS function for identifying spatial overlay:
+
+```{r}
+sel_africa = st_covered_by(world_centroids, africa, sparse = FALSE)
+summary(sel_africa)
+```
+
+This shows that there are 56 countries in Africa.
+We can check if they are the same countries as follows:
+
+```{r}
+africa_centroids2 = world_centroids[sel_africa,]
+identical(africa_centroids, africa_centroids2)
+```
+
+## Exercises
+
+- Perform the same operations and map making for another continent of your choice.
+- Bonus: Download some global geographic data and add attribute variables assigning them to the continents of the world.
+
+
 # References
diff --git a/code/01-introduction.R b/code/01-introduction.R
index c8ac0cc23..83aab3d8f 100644
--- a/code/01-introduction.R
+++ b/code/01-introduction.R
@@ -15,8 +15,8 @@ world = st_read(f)
 ## ------------------------------------------------------------------------
 class(world)
 
-## ------------------------------------------------------------------------
-world[1:2, 1:3]
+## ---- results='hide'-----------------------------------------------------
+summary(world)
 
 ## ------------------------------------------------------------------------
 world_sp = as(object = world, Class = "Spatial")
diff --git a/code/02-read-write.R b/code/02-read-write.R
index cd46b341e..837d43589 100644
--- a/code/02-read-write.R
+++ b/code/02-read-write.R
@@ -1,3 +1,7 @@
+## ---- echo=FALSE, include=FALSE------------------------------------------
+if(!exists("world"))
+        source("code/01-introduction.R")
+
 ## ---- results='hide'-----------------------------------------------------
 library(microbenchmark)
 bench_read = microbenchmark(times = 5,
@@ -9,23 +13,24 @@ bench_read = microbenchmark(times = 5,
 bench_read$time[1] / bench_read$time[2]
 
 ## ---- echo=FALSE, results='hide'-----------------------------------------
-w_files = list.files(pattern = "w\\.")
-file.remove(w_files)
+world_files = list.files(pattern = "world\\.")
+file.remove(world_files)
 
 ## ---- warning=FALSE, results='hide'--------------------------------------
 bench_write = microbenchmark(times = 1,
-        st_write(w, "w.geojson"),
-        st_write(w, "w.shp"),
-        st_write(w, "w.gpkg")
+        st_write(world, "world.geojson"),
+        st_write(world, "world.shp"),
+        st_write(world, "world.gpkg")
 )
 
 ## ---- echo=FALSE, results='hide'-----------------------------------------
-w_files = list.files(pattern = "w\\.")
-file.remove(w_files)
+world_files = list.files(pattern = "world\\.")
+file.remove(world_files)
 
 ## ------------------------------------------------------------------------
 bench_write
 
 ## ------------------------------------------------------------------------
-st_drivers()[1:2,]
+sf_drivers = st_drivers()
+head(sf_drivers, n = 2)
 
diff --git a/code/03-attr.R b/code/03-attr.R
index 4afd3ec23..06668e2ba 100644
--- a/code/03-attr.R
+++ b/code/03-attr.R
@@ -1,5 +1,6 @@
 ## ---- echo=FALSE, include=FALSE------------------------------------------
-source("code/01-introduction.R")
+if(!exists("world"))
+        source("code/01-introduction.R")
 
 ## ------------------------------------------------------------------------
 class(world)
@@ -15,60 +16,59 @@ st_geometry(world_df) = NULL
 class(world_df)
 
 ## ---- eval=FALSE---------------------------------------------------------
-#> world[1:6, ] # subsetting rows
+#> world[1:6, ] # subset rows
 
-## ------------------------------------------------------------------------
-world[, 1:3] # subsetting columns
+## ---- eval=FALSE---------------------------------------------------------
+#> world[, 1:3] # subset columns
 
 ## ---- message=FALSE------------------------------------------------------
 library(dplyr)
 
 ## ------------------------------------------------------------------------
-world = select(world, name, continent, population = pop_est)
-head(world)
+world_orig = world # create copy of world dataset for future reference
+world = select(world_orig, name, continent, population = pop_est)
+head(world, n = 2)
 
 ## ---- eval=FALSE---------------------------------------------------------
-#> world = world[c("name", "continent", "pop_est")] # subset columns by name
-#> names(world)[3] = "population" # rename column manually
+#> world2 = world_orig[c("name", "continent", "pop_est")] # subset columns by name
+#> names(world2)[3] = "population" # rename column manually
 
 ## ------------------------------------------------------------------------
-world_few_cols2 = world %>%
-        select(., name, continent)
-
-head(world_few_cols2)
+world4 = world_orig %>%
+        select(name, continent)
 
 ## ------------------------------------------------------------------------
 # ==, !=, >, >=, <, <=, &, |
 
 # subsetting simple feature rows by values
-world_few_rows = world[world$pop_est>1000000000, ]
+world_few_rows = world[world$population > 1000000000,]
 
 #OR
 world_few_rows = world %>% 
-        filter(., pop_est>1000000000)
+        filter(population > 1000000000)
 
 head(world_few_rows)
 
 ## ------------------------------------------------------------------------
 # add a new column
 world$area = raster::area(as(world, "Spatial")) / 1000000 #it there any function for area calculation on sf object?
-world$pop_density = world$pop_est / world$area
+world$pop_density = world$population / world$area
 
 # OR
 
 world = world %>% 
-        mutate(area=raster::area(as(., "Spatial")) / 1000000) %>% 
-        mutate(pop_density=pop_est/area)
+        mutate(area = raster::area(as(., "Spatial")) / 1000000) %>%
+        mutate(pop_density = population / area)
 
-## ------------------------------------------------------------------------
-# data summary
+## ---- results='hide'-----------------------------------------------------
+# data summary (not shown)
 summary(world)
 
-# data summary by groups
+# data summary by groups (not shown)
 world_continents = world %>% 
         group_by(continent) %>% 
-        summarise(continent_pop=sum(pop_est), country_n=n())
-world_continents
+        summarise(continent_pop = sum(population), country_n = n())
+                  world_continents
 
 ## ------------------------------------------------------------------------
 # sort variables
@@ -87,6 +87,6 @@ class(world_st)
 # OR
 
 world_st2 = world
-world_st2 = world_st2 %>% st_set_geometry(., NULL)
+world_st2 = world_st2 %>% st_set_geometry(NULL)
 class(world_st2)
 
diff --git a/code/04-basic-map.R b/code/04-basic-map.R
index e2628f148..80c6d95c1 100644
--- a/code/04-basic-map.R
+++ b/code/04-basic-map.R
@@ -1,4 +1,38 @@
+## ---- echo=FALSE, include=FALSE------------------------------------------
+if(!exists("world"))
+        source("code/03-attr.R")
+
 ## ----sfplot, fig.cap="Plotting with sf, with multiple variables (left) and a single variable (right).", out.width="49%", fig.show='hold', warning=FALSE----
-plot(w)
-plot(w["pop_est"])
+plot(world)
+plot(world["population"])
+
+## ---- fig.show='hide'----------------------------------------------------
+plot(world["pop_density"])
+africa = world_continents[1,]
+plot(africa, add = TRUE, col = "red")
+
+## ----nigeria, warning=FALSE, echo=FALSE, fig.cap="Map of Nigeria in context illustrating sf's plotting capabilities"----
+nigeria = filter(world, name == "Nigeria")
+bb_africa = st_bbox(africa)
+plot(africa[2], col = "white", lwd = 3, main = "Nigeria in context", border = "lightgrey")
+# plot(world, lty = 3, add = TRUE, border = "grey")
+plot(world, add = TRUE, border = "grey")
+plot(nigeria, col = "yellow", add = TRUE, border = "darkgrey")
+ncentre = st_centroid(nigeria)
+ncentre_num = st_coordinates(ncentre)
+text(x = ncentre_num[1], y = ncentre_num[2], labels = "Nigeria")
+
+## ---- out.width="50%", fig.cap="Centroids in Africa"---------------------
+world_centroids = st_centroid(world)
+plot(world_centroids[1])
+africa_centroids = world_centroids[africa,]
+plot(africa_centroids, add = TRUE, cex = 2)
+
+## ------------------------------------------------------------------------
+sel_africa = st_covered_by(world_centroids, africa, sparse = FALSE)
+summary(sel_africa)
+
+## ------------------------------------------------------------------------
+africa_centroids2 = world_centroids[sel_africa,]
+identical(africa_centroids, africa_centroids2)
 
diff --git a/index.Rmd b/index.Rmd
index 9b05e241f..aaac30dbc 100644
--- a/index.Rmd
+++ b/index.Rmd
@@ -21,7 +21,6 @@ url: https\://Robinlovelace.github.io/rgeobook/
 # source("code/before_script.R")
 ```
 
-
 # Prerequisites {-}
 
 To run the following code you must have **sf** installed, e.g. via: