Chapter 10 A Tidy Approach to Spatial Data: Simple Features

Note that if you, the student, wish to run all the code yourself, you should download the inputs directory as a zipped file by going here with a web browser and then clicking the big “Download” button on the right. Once you have downloaded that, unzip it and put the whole inputs directory in the current working directory where you are working with R.

10.1 Prelims

If you don’t already have it:

install.packages("sf")

Also, as of June 5, 2017, this requires the development version of ggplot2:

devtools::install_github("tidyverse/ggplot2")

Be warned that this seems to break some of the earlier ggmap.

library(tidyverse)
library(sf)
## Linking to GEOS 3.4.2, GDAL 2.1.2, proj.4 4.9.1
library(raster)

10.2 Reading spatial data into sf objects

Let’s just quickly put the code in to reproduce most of the stuff from the previous lecture, but using sf.

Note that there are two functions to read shapefiles: st_read() and read_sf(). The latter seems preferable to me as it does not read stringsAsFactors ever.

streams <- read_sf("inputs/santa-cruz-county/Streams/Streams.shp")
watersheds <- read_sf("inputs/santa-cruz-county/Watersheds/Watersheds.shp")
lights <- read_sf("inputs/santa-cruz-county/Street_Lights/Street_Lights.shp")

Have a look at one of them. Note that the whole thing is a data frame, so you can directly as_tibble() it:

as_tibble(streams)
## # A tibble: 1,514 x 9
##    OBJECTID       STREAM_NM   STREAMTYPE  GNIS_ID    ID ENABLED Named
##       <dbl>           <chr>        <chr>    <chr> <chr>   <chr> <chr>
##  1        1     Adams Creek INTERMITTENT     <NA>     0    <NA>   Yes
##  2        2      Alba Creek    PERENNIAL 00218063     1    <NA>   Yes
##  3        3     Amaya Creek    PERENNIAL 00218217     2    <NA>   Yes
##  4        4 Ano Nuevo Creek    PERENNIAL 00254567     3    <NA>   Yes
##  5        5     Aptos Creek    PERENNIAL 00254571     4    <NA>   Yes
##  6        6 Archibald Creek INTERMITTENT 00218350     5    <NA>   Yes
##  7        7   Baldwin Creek    PERENNIAL 00218614     6    <NA>   Yes
##  8        8     Bates Creek    PERENNIAL 00218726     7    <NA>   Yes
##  9        9      Bean Creek    PERENNIAL 00218782     8    <NA>   Yes
## 10       10   Bennett Creek    PERENNIAL 00219040     9    <NA>   Yes
## # ... with 1,504 more rows, and 2 more variables: SHAPESTLen <dbl>,
## #   geometry <simple_feature>

10.3 Plotting sf objects

In ggplot this is made easy with the geom_sf() which knows what type of simple feature you have, and plots it appropriately.

You also add coord_sf() which will deal with projections, on the fly, and also takes care of aspect ratios and draws a graticule.

ggplot() +
  geom_sf(mapping = aes(colour = STREAMTYPE), data = streams) +
  coord_sf()

Check it out, there is a graticule. Note that the long-lat lines are not perfectly square. That is because this object is not in long-lat:

st_crs(streams)
## $epsg
## [1] NA
## 
## $proj4string
## [1] "+proj=lcc +lat_1=37.06666666666667 +lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000 +y_0=500000.0000000002 +datum=NAD83 +units=us-ft +no_defs"
## 
## attr(,"class")
## [1] "crs"

Aha! See, it is in lcc = lambert conic conformal.

What happens if we try to put an open street maps on this? It fails…

How about if we try to use spraster_rgb to put a raster background in there? That is a little complicated. So, instead let’s have a look at interactive mapping later.

But before that I want to briefly touch on geometrical intersections.

10.4 Geometric operations

Just briefly, remember, we wondered last week how to find all the streams that are in the Scott Creek watershed. This can be done with the raster package, nicely. But with the sf package we can also use st_intersection() to spatially intersect streams with watersheds. When we do this, the cool thing is that all the columns of data associated with each get joined just as you would hope they do!!

streamsheds <- st_intersection(streams, watersheds)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

Those columns get joined on the basis of the stream being in the watershed. Cool!

Check it out:

as_tibble(streamsheds)
## # A tibble: 1,522 x 15
##    OBJECTID         STREAM_NM   STREAMTYPE  GNIS_ID    ID ENABLED Named
##  *    <dbl>             <chr>        <chr>    <chr> <chr>   <chr> <chr>
##  1        4   Ano Nuevo Creek    PERENNIAL 00254567     3    <NA>   Yes
##  2       25     Cascade Creek    PERENNIAL 00220660    24    <NA>   Yes
##  3       38      Elliot Creek    PERENNIAL 00223153    37    <NA>   Yes
##  4       41      Finney Creek    PERENNIAL 00223512    40    <NA>   Yes
##  5       47       Gazos Creek    PERENNIAL     <NA>    46    <NA>   Yes
##  6       51  Green Oaks Creek    PERENNIAL 00224554    50    <NA>   Yes
##  7       91  Old Womans Creek    PERENNIAL 00229991    90    <NA>   Yes
##  8      129 White House Creek    PERENNIAL     <NA>   128    <NA>   Yes
##  9      131      Willow Gulch INTERMITTENT     <NA>   130    <NA>   Yes
## 10      162       Stream 1141        SWALE     <NA>   161    <NA>    No
## # ... with 1,512 more rows, and 8 more variables: SHAPESTLen <dbl>,
## #   OBJECTID.1 <dbl>, NAME <chr>, COUNT_ <dbl>, MapKey <chr>,
## #   SHAPESTAre <dbl>, SHAPESTLen.1 <dbl>, geometry <simple_feature>

Plot it with different colors for different watersheds to see what I mean:

ggplot(streamsheds) +
  geom_sf(aes(colour = NAME)) +
  coord_sf()

OK, now let’s make some interactive maps.

10.5 Interactive plotting with leaflet via the mapview package

This is truly mind-blowing. The authors of the mapview package tell us:

mapview is an R package created to help researchers during their spatial data analysis workflow. It provides functions to very quickly and conveniently create interactive visualisations of spatial data. It was created to fill the gap of quick (not presentation grade) interactive plotting to examine and visually investigate both aspects of spatial data, the geometries and their attributes.”

As such, it appears that customization of certain things might not be easily-accessible, but it is so easy and still so customizable that it is really, really amazing.

Get it from CRAN:

install.packages("mapview")

and then be sure to load it:

library(mapview)
## Loading required package: leaflet

The main function is mapView which does multiple dispatch on a huge variety of different object types (rasters, points, polygons, etc), including sf classes. It creates an interactive web-map powered by leaflet, an open-source JavaScript library for mobile-friendly interactive maps.

Let’s look at it in action:

mapView(streams)