Skip to content

Commit

Permalink
Many updates - ready for teaching
Browse files Browse the repository at this point in the history
  • Loading branch information
Robinlovelace committed Mar 28, 2017
1 parent 0b1722f commit bc2e0e1
Show file tree
Hide file tree
Showing 8 changed files with 150 additions and 42 deletions.
6 changes: 3 additions & 3 deletions 01-introduction.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
6 changes: 3 additions & 3 deletions 03-attr.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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())
Expand Down
70 changes: 70 additions & 0 deletions 04-basic-map.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 2 additions & 2 deletions code/01-introduction.R
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
21 changes: 13 additions & 8 deletions code/02-read-write.R
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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)

46 changes: 23 additions & 23 deletions code/03-attr.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
## ---- echo=FALSE, include=FALSE------------------------------------------
source("code/01-introduction.R")
if(!exists("world"))
source("code/01-introduction.R")

## ------------------------------------------------------------------------
class(world)
Expand All @@ -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
Expand All @@ -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)

38 changes: 36 additions & 2 deletions code/04-basic-map.R
Original file line number Diff line number Diff line change
@@ -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)

1 change: 0 additions & 1 deletion index.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down

0 comments on commit bc2e0e1

Please sign in to comment.