30

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

`

library(sf)

## Linking to GEOS 3.5.1, GDAL 2.1.3, PROJ 4.9.2

library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

## ✔ ggplot2 3.1.0       ✔ purrr   0.3.2  
## ✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0

## ── Conflicts ───────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

library(spData)

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

summary(st_geometry(world))

##  MULTIPOLYGON     epsg:4326 +proj=long... 
##           177             0             0

#>  MULTIPOLYGON     epsg:4326 +proj=long... 
#>           177             0             0

old_par = par() # save graphical parameters

# generate version based on code in the book
plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)

## Warning in st_centroid.sf(world, of_largest = TRUE): st_centroid assumes
## attributes are constant over geometries of x

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

plot(st_geometry(world_cents), add = TRUE, cex = cex, reset = TRUE)

# generate version in the book from online script:
par(old_par)

## Warning in par(old_par): graphical parameter "cin" cannot be set

## Warning in par(old_par): graphical parameter "cra" cannot be set

## Warning in par(old_par): graphical parameter "csi" cannot be set

## Warning in par(old_par): graphical parameter "cxy" cannot be set

## Warning in par(old_par): graphical parameter "din" cannot be set

## Warning in par(old_par): graphical parameter "page" cannot be set

s = "https://github.com/Robinlovelace/geocompr/raw/master/code/02-contpop.R"
source(s)

## Warning in st_centroid.sf(world_proj, of_largest_polygon = TRUE):
## st_centroid assumes attributes are constant over geometries of x

Morocco = world[world$name_long == "Morocco", ]
africa = world[world$continent == "Africa", ]
plot(st_geometry(Morocco), col = "white", lwd = 3, main = "Morocco in context", border = "lightgrey", expandBB = c(0.5, 0.2, 0.5, 0.2))
plot(st_geometry(world), lty = 3, add = TRUE, border = "grey")
plot(st_geometry(Morocco), col = "red", add = TRUE, border = "darkgrey")
a = africa[grepl("Morocco", africa$name_long), ]
ncentre = st_centroid(a)

## Warning in st_centroid.sf(a): st_centroid assumes attributes are constant
## over geometries of x

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

ncentre_num = st_coordinates(ncentre)
text(x = ncentre_num[, 1], y = ncentre_num[, 2], labels = a$name_long)

set.seed(2018-09-03)

library(sf)
library(raster)
library(tidyverse)
library(spData)
#library(spDataLarge)

world_sfc = st_geometry(world)
world_sfc_mirror = world_sfc * c(1, -1)

## Warning in mapply(function(x, y) {: longer argument not a multiple of
## length of shorter

#> Warning in mapply(function(x, y) {: longer argument not a multiple of
#> length of shorter
plot(world_sfc)

library(sf)
library(raster)
library(tidyverse)
library(spData)
library(tmap)
library(leaflet)
library(cartogram)

africa = world %>% 
  filter(continent == "Africa", !is.na(iso_a2)) %>% 
  left_join(worldbank_df, by = "iso_a2") %>% 
  dplyr::select(name, subregion, gdpPercap, HDI, pop_growth) %>% 
  st_transform("+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25")

plot(africa["HDI"])

tm_shape(africa) + tm_polygons("HDI")

library(ggplot2)
ggplot(africa) +
  geom_sf(aes(geometry = geom, fill = HDI)) 

library(leaflet)
leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = st_transform(africa, 4326))

library(cartogram)
africa_carto = cartogram_cont(africa, "HDI", itermax = 1)

## Warning in cartogram_cont.sf(africa, "HDI", itermax = 1): NA not allowed in
## weight vector. Features will be removed from Shape.

## Mean size error for iteration 1: 6.03537025368655

#> Warning in cartogram_cont.sf(africa, "HDI", itermax = 5): NA not allowed in
#> weight vector. Features will be removed from Shape.
#> Mean size error for iteration 1: 6.03537025368671
#> Mean size error for iteration 2: 3.64410634493783
#> Mean size error for iteration 3: 2.45300713828664
#> Mean size error for iteration 4: 1.85135460746932
#> Mean size error for iteration 5: 1.53025529469271
tm_shape(africa_carto) + tm_polygons("HDI")

# >0.8 - Very high human development
# >0.7 - High human development
# >0.55 - Medium human development
# <0.55 - Low human development
sm1 = tm_shape(africa) +
  tm_polygons(
  col = "HDI",
  title = "Human Development Index",
  breaks = c(0, 0.55, 0.7, 0.8),
  labels = c("Low", "Medium", "High"),
  palette = "YlGn")
sm1

sm2 = tm_shape(africa) +
  tm_polygons(col = "subregion",
              title = "Subregion", 
              palette = "Set2")
tmap_arrange(sm2, sm1)

northern_africa = filter(africa, subregion == "Northern Africa")
tm_shape(northern_africa) +
  tm_polygons(col = c("HDI", "pop_growth")) +
  qtm(africa, fill = NULL)

tm_shape(northern_africa) +
  tm_polygons("pop_growth", style = "jenks", palette = "BuPu") +
  tm_facets(by = "name", drop.NA.facets = TRUE)

# with tmap:
tmap_mode("view")

## tmap mode set to interactive viewing

#> tmap mode set to interactive viewing
tm_shape(northern_africa) +
  tm_polygons("pop_growth") +
  tm_scale_bar() # bonus scale bar

tmap_mode("plot")

## tmap mode set to plotting

#> tmap mode set to plotting

# with mapview (auto scale + legend)
library(mapview)
mapview(northern_africa, zcol = "pop_growth")