dgl5gw/geocompr
Geocomputation with R: a book

Package index
Search the dgl5gw/geocompr package
Vignettes
Functions
0
Source code
34
Man pages
0
Home
/
GitHub
/
dgl5gw/geocompr
/

In dgl5gw/geocompr: Geocomputation with R: a book 

options(htmltools.dir.version = FALSE)
knitr::opts_chunk$set(fig.align = "center", warning = FALSE, message = FALSE)

Prerequisites

install.packages(c("raster", "sf", "rasterVis", "tmap", "leaflet",
                   "widgetframe", "tidyverse", "spData", "geofacet",
                   "tidycensus", "linemap", "gapminder", "tigris", 
                   "osmdata", "mapview", "opencage"))

Prerequisites

Note: tidycensus and opencage require API keys.

To get a tidycensus API key

Go to http://api.census.gov/data/key_signup.html and complete the form there. Replace text YOUR API KEY GOES HERE with your census API key. Execute the code below:

library(tidycensus)
census_api_key("YOUR API KEY GOES HERE", install = TRUE)

To get a opencage API key

Go to https://geocoder.opencagedata.com/pricing, select the Free Trial, and complete the form there. Replace text YOUR API KEY GOES HERE with your opencage API key. Execute the code below:

cat("OPENCAGE_KEY=YOUR API CODE HERE", file=file.path(normalizePath("~/"), ".Renviron"), append=TRUE)

Prerequisites

Check if the keys are installed

Restart R/RStudio (Session -> Restart R) and execute the code below:

Sys.getenv("CENSUS_API_KEY")
Sys.getenv("OPENCAGE_KEY")

You should see your keys print to screen.

More info can be found here and here.

class: center, middle

Part I - Examples

Static visualization (I)

knitr::include_graphics("figs/01_world_globe.png")

Static visualization (II)

knitr::include_graphics("figs/02_state_unemp.png")

Static visualization (III)

knitr::include_graphics("figs/03_ohio_linemap.png")

Animation

knitr::include_graphics("figs/04_anim.gif")

Interactive visualization

Remote sensing calculations

knitr::include_graphics("figs/06_ndvi.png")

Terrain calculations

knitr::include_graphics("figs/07_terrain.png")

Geocoding

What's the coordinates of New York, Los Angeles, Chicago, Houston, and Philadelphia?

class: center, middle

Part II - An intro to R

How to get R

Rstudio

How to get help?

# if you know a function name
?crop

Online help

How to start with R?

R packages

install.packages("raster")

library(raster)

class: center, middle

Part III - An intro to spatial R

Spatial R

d = readr::read_csv("figs/gis-vs-gds-table.csv")
knitr::kable(x = d, caption = "Differences in emphasis between the fields of Geographic Information Systems (GIS) and Geographic Data Science (GDS)", format = "html")

The sf package

The sf package in an R implementation of Simple Features. This package incorporates: - A new spatial data class system in R - Functions for reading and writing data - Tools for spatial operations on vectors

Most of the functions in this package starts with prefix st_.

devtools::install_github("edzer/sfr") # development version

or

install.packages("sf") # stable version

You need a recent version of the GDAL, GEOS, Proj.4, and UDUNITS libraries installed for this to work on Mac and Linux. More information on that at https://github.com/edzer/sfr.

library(sf)

Reading spatial data

wrld = st_read("data/wrld.shp")

ham = st_read("data/hamilton_county.gpkg")

Reading spatial data - text files

ham_cities_sf = st_read("data/hamiltion_cities.csv", 
                        options = c("X_POSSIBLE_NAMES=X",
                                    "Y_POSSIBLE_NAMES=Y",
                                    "KEEP_GEOM_COLUMNS=NO"))

Writing spatial data

file.remove(c("data/new_wrld.shp", "data/new_wrld.gpkg"))

st_write(wrld, "data/new_wrld.shp")

st_write(wrld, "data/new_wrld.gpkg")

file.remove(c("data/new_wrld.shp", "data/new_wrld.gpkg"))

sf structure

The sf object is a table (data.frame) with an additional geometry column - typically named geom or geometry and spatial metadata (geometry type, dimension, bbox, epsg (SRID), proj4string).

wrld[1:5, 1:3]

Attributes - the dplyr package

library(dplyr)

wrld_fil = filter(wrld, pop < 297517)

wrld_fil[1:3, ]

Attributes - the dplyr package

wrld_mut = mutate(wrld, pop_density = pop/area_km2)

wrld_mut[1:3, ]

Attributes - the dplyr package

wrld_sum1 = summarize(wrld, pop_sum = sum(pop, na.rm = TRUE), 
                      pop_mean = mean(pop, na.rm = TRUE), 
                      pop_median = median(pop, na.rm = TRUE))

wrld_sum1[1, ]

Attributes - the dplyr package

wrld_sum1 = wrld %>% 
        group_by(continent) %>% 
        summarize(pop_sum = sum(pop, na.rm = TRUE), 
                      pop_mean = mean(pop, na.rm = TRUE), 
                      pop_median = median(pop, na.rm = TRUE))

wrld_sum1[1:3, ]

CRS and Reprojection

st_crs(wrld)

wrld_3857 = st_transform(wrld, 3857)
st_crs(wrld_3857)

CRS and Reprojection

png("figs/coord_compare.png", width = 800, height = 500)
par(mfrow = c(1, 2), mar=c(0,0,0,0))
plot(wrld[0]);plot(wrld_3857[0])
dev.off()

Spatial operations

continents = wrld %>% 
        st_set_precision(1000) %>% 
        st_union()
plot(continents)

continents = wrld %>% 
        st_set_precision(1000) %>% 
        st_union()
par(mar = c(0, 0, 0, 0))
plot(continents)

The raster package

The raster package consists of method and classes for raster processing. It allows to:

Basic raster operations are illustrated at https://rpubs.com/etiennebr/visualraster.

This package has three object classes:

library(raster)

Reading

dem = raster("data/srtm.tif")
dem

Writing

file.remove(c("data/new_dem.tif", "data/new_dem2.tif"))

writeRaster(dem, "data/new_dem.tif")

writeRaster(dem, "data/new_dem2.tif", 
            datatype = "FLT4S", options=c("COMPRESS=DEFLATE"))

file.remove(c("data/new_dem.tif", "data/new_dem2.tif"))

writeFormats()

raster structure

dem

inMemory(dem)

Attributes

values_dem = getValues(dem)

values_dem[1:50]

Attributes

.pull-left[

new_dem = dem + 50

plot(new_dem)

]

.pull-right[

new_dem2 = dem * new_dem

plot(new_dem2)

]

CRS and Reprojection

crs(dem)

dem3857 = projectRaster(dem, 
                         crs="+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0
                         +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0
                         +units=m +no_defs")

http://spatialreference.org/

CRS and Reprojection

png("figs/coord_compare_raster.png", width = 800, height = 400)
par(mfrow = c(1, 2), mar=c(4,4,4,8))
# par(mfrow = c(1, 2))
plot(dem);plot(dem3857)
dev.off()

Extract

ham_cities_sp = as(ham_cities_sf, "Spatial")
raster::extract(dem, ham_cities_sp)

ham_cities_sf$dem = raster::extract(dem, ham_cities_sp)
ham_cities_sf

Crop

ham84 = st_transform(ham, 4326)
ham_sp = as(ham84, "Spatial")

.pull-left[ 

dem_crop = crop(dem, ham_sp)

plot(dem_crop)

]

.pull-right[ 

dem_mask = mask(dem_crop, ham_sp)

plot(dem_mask)

]

Reading spatial data - data packages

datapackages = tibble::tribble(~`Package name`, ~Description,
                               "osmdata", "Download and import of OpenStreetMap data",
                               "raster", "The getData() function downloads and imports administrative country, SRTM/ASTER elevation, WorldClim data",
                               "rnaturalearth", "Functions to download Natural Earth vector and raster data, including world country borders",
                               "tigirs", "Census TIGER/Line shapefiles in R",
                               "USAboundaries", "Contemporary state, county, and congressional district boundaries for the United States of America, as well as historical boundaries from 1629 to 2000 for states and counties",
                               "rnoaa", "An interface to many NOAA data sources",
                               "rWBclimate", "The World Bank climate data")
knitr::kable(datapackages, format = "html")

Basic maps

plot(wrld[0])

plot(wrld["pop"])

png("figs/plot_compare.png", width = 800, height = 300)
par(mfrow = c(1, 2), mar=c(0,0,1,0))
plot(wrld[0]);plot(wrld["pop"])
dev.off()

Basic maps

plot(dem)

rasterVis

library(rasterVis)
my_theme = rasterTheme(region=brewer.pal("RdYlGn", n = 9))
p = levelplot(dem_crop, margin = FALSE, par.settings = my_theme)
p = p + layer(sp.lines(ham_sp, lwd = 3, col = "darkgrey"))
p + layer(sp.points(ham_cities_sp, pch = 19, col = "black"))

tmap

library(tmap)
tm_shape(wrld, projection="wintri") +
        tm_polygons("lifeExp", style="pretty", palette="RdYlGn",
                    auto.palette.mapping=FALSE, title=c("Life expectancy")) +
        tm_style_grey()

leaflet

library(leaflet)
leaflet(ham_sp) %>%
        addPolygons() %>% 
        addMarkers(data=ham_cities_sp, popup=~as.character(name))

library(widgetframe)
library(leaflet)
l = leaflet(ham_sp) %>%
        addPolygons() %>% 
        addMarkers(data=ham_cities_sp, popup=~as.character(name))
frameWidget(l, height = "400")

class: center, middle

Part IV - Examples

Static visualization (I)

knitr::include_graphics("figs/01_world_globe.png")

Static visualization (I)

library(tidyverse)
library(sf)
library(spData)

world = st_read(system.file("shapes/world.gpkg", package = "spData"))
world_globe = sf::st_transform(world,
        "+proj=ortho +x_0=0 +y_0=0 +lon_0=-100 +lat_0=40")

p1 = ggplot(world_globe) +
        geom_sf(data = world_globe, aes(fill = lifeExp))  +
        scale_fill_viridis_c(name = "Life expectancy: ") +
        theme_void()
p1
ggsave(p1, file="figs/01_world_globe.png", 
       width = 70, height = 60, units = "mm", scale = 2)

Static visualization (II)

knitr::include_graphics("figs/02_state_unemp.png")

Static visualization (II)

https://hafen.github.io/geofacet/

library(geofacet)
library(ggplot2)
data("state_unemp", package = "geofacet")
p2 = ggplot(state_unemp, aes(year, rate)) +
        geom_line() +
        facet_geo(~state, grid = "us_state_grid2", label = "name") +
        scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) +
        labs(title = "Seasonally Adjusted US Unemployment Rate 2000-2016",
             caption = "Data Source: bls.gov",
             x = "Year",
             y = "Unemployment Rate (%)") +
        theme(strip.text.x = element_text(size = 2)) + 
        theme_bw()
ggsave(p2, file="figs/02_state_unemp.png", 
       width = 70, height = 60, units = "mm", scale = 4)

Static visualization (III)

knitr::include_graphics("figs/03_ohio_linemap.png")

Static visualization (III)

https://github.com/rCarto/linemap

library(sf)
library(tidycensus)
library(linemap)
# A key can be acquired at http://api.census.gov/data/key_signup.html
# census_api_key("your key")
oh_pop = get_acs(geography = "tract", variables = "B01003_001E", 
                  state = "OH", geometry = TRUE) 
plot(oh_pop$geometry)
oh_pop
oh_pop2 = st_transform(oh_pop, 26917)

oh_pop2_grid = getgrid(x = oh_pop2, cellsize = 2500, var = "estimate")

png("figs/03_ohio_linemap.png", width = 650, height = 600)
opar = par(mar = c(0,0,0,0))
plot(st_geometry(oh_pop2), col="lightsteelblue3", border = NA, bg = "lightsteelblue1")
linemap(x = oh_pop2_grid, var = "estimate", k = 3, threshold = 1,
        col = "lightsteelblue3", border = "white", lwd = 0.8,
        add = TRUE)
par(opar)
dev.off()

Animation

knitr::include_graphics("figs/04_anim.gif")

Animation

library(gapminder)
library(tidyverse)
library(sf)
library(tmap)
library(spData)
europe_area = st_polygon(list(rbind(c(-13, 34), c(-13, 72),
                                    c(43, 72), c(43, 34), c(-13, 34)))) %>% 
        st_sfc(crs = st_crs(world)) 
europe = st_intersection(world, europe_area) 
europe_gdp = europe %>% 
        inner_join(gapminder, by = c("name_long" = "country"))
m1 = tm_shape(europe) + 
        tm_fill("grey") + 
     tm_shape(europe_gdp) + 
        tm_fill("gdpPercap.y", title = "GDP per capita: ") +
        tm_borders() + 
        tm_facets(by = "year", nrow = 1, ncol = 1, drop.units = TRUE) + 
        tm_legend(text.size = 1, title.size = 1.2, 
                  position = c("left", "TOP"), height = 0.3)

animation_tmap(m1, filename = "figs/04_anim.gif", width = 650, height = 600)

library(gapminder)
library(tidyverse)
library(sf)
library(tmap)
library(spData)
world_gdp = world %>% 
        left_join(gapminder, by = c("name_long" = "country"))

m1 = tm_shape(world_gdp, bbox = raster::extent(-13, 43, 34, 72)) + 
        tm_polygons("grey") +
     tm_shape(world_gdp) + 
        tm_polygons("gdpPercap.y", title = "GDP per capita: ", breaks = seq(0, 25000, by = 5000)) +
        tm_facets(by = "year", nrow = 1, ncol = 1, drop.units = TRUE, showNA = FALSE) + 
        tm_legend(text.size = 1,
                  title.size = 1.2,
                  position = c("left", "TOP"),
                  height = 0.3)

animation_tmap(m1, filename = "figs/04_anim.gif", width = 650, height = 600)

Interactive visualization

Interactive visualization

library(tigris)
library(leaflet)
library(sf)
library(tidyverse)
library(osmdata)
library(widgetframe)
oh = counties(state="OH")
hamilton = oh %>% 
        st_as_sf() %>% 
        filter(NAME == "Hamilton") %>% 
        st_transform(4326)
pubs = opq(bbox = "Cincinnati OH") %>% 
  add_osm_feature(key = "amenity", value = "pub") %>% 
  osmdata_sf()
l1 = leaflet(hamilton) %>%
        addTiles(group = "OSM") %>%
        addProviderTiles(providers$Stamen.Watercolor, group = "Watercolor") %>% 
        addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>% 
        addLayersControl(    baseGroups = c("OSM", "Watercolor", "Toner Lite")) %>% 
        addPolygons() %>% 
        addMarkers(data = pubs$osm_points, popup=~as.character(name))
saveWidget(l1, file = "05_leaflet.html", selfcontained = FALSE)

Remote sensing calculations

knitr::include_graphics("figs/06_ndvi.png")

Remote sensing calculations

library(raster)
# download_landsat8 = function(destination_filename, band){
#         filename = paste0("http://landsat-pds.s3.amazonaws.com/L8/038/034/LC80380342015230LGN00/LC80380342015230LGN00_B", band, ".TIF")
#         download.file(filename, destination_filename, method='auto')
# }
# download_landsat8("data/landsat_b4.tif", band = 4)
# download_landsat8("data/landsat_b5.tif", band = 5)
b4 = raster("data/landsat_b4.tif")
b5 = raster("data/landsat_b5.tif")

ndvi_calc = function(b4, b5){
        ndvi = (b5 - b4)/(b5 + b4)
        ndvi
}

my_ndvi = overlay(b4, b5, fun = ndvi_calc)

png("figs/06_ndvi.png", width = 650, height = 500)
par(mfrow=c(1,2))
plot(my_ndvi, main = "NDVI")
hist(my_ndvi, main = "", xlab = "NDVI")
dev.off()

Terrain calculations

knitr::include_graphics("figs/07_terrain.png")

Terrain calculations

library(raster)
# getData('SRTM', lon = 5, lat = 45, path = "data")
srtm = raster("data/srtm_38_04.tif")
my_area = extent(c(8, 10, 41.3, 43.1))
my_srtm = crop(srtm, my_area)
plot(my_srtm)

library(mapview)
mapView(my_srtm)

my_terrain = terrain(my_srtm, opt = c("slope", "aspect", "tpi", "tri", "roughness", "flowdir"))

png("figs/07_terrain.png", width = 650, height = 500)
plot(my_terrain)
dev.off()

Geocoding

What's the coordinates of New York, Los Angeles, Chicago, Houston, and Philadelphia?

Geocoding

# http://happygitwithr.com/api-tokens.html
library(opencage)
library(tidyverse)
library(sf)
library(mapview)
library(htmlwidgets)

my_places = c("New York", "Los Angeles", "Chicago",
              "Houston", "Philadelphia")

output = my_places %>%
        map(opencage_forward, limit = 1) %>% 
        map_df(1)

output_sf = st_as_sf(output, 
                     coords = c("geometry.lng", "geometry.lat"), crs = 4326)
plot(output_sf$geometry, axes = TRUE)

m1 = mapview(output_sf)
mapshot(m1, "08_geocoding.html")

class: center, middle

Part V - The end

The online version of the book is at http://robinlovelace.net/geocompr/ and its source code at https://github.com/robinlovelace/geocompr.

We encourage contributions on any part of the book, including:

Please see our_style.md for the books style.

References

References

References



dgl5gw/geocompr documentation built on May 18, 2019, 8:11 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close