GIS in R

Jared P. Lander

• Let’s talk about maps
• Ok, this is a globe
• Flatten into a map

All About Relations

• All About Relations
• Relations between what?

• Points
• Lines
• Polygons

• Points
• Lines
• Polygons
• Multi-versions
• How they relate to each other

Today’s Goals

• What we will do today

• Draw Maps with Different Packages
• Compute Relations Between Points and Polygons
• Solve a Fun Mapping Puzzle

• Draw Maps with Different Packages
• Compute Relations Between Points and Polygons
• Solve a Fun Mapping Puzzle
• Overlap between topics

Maps

• Different ways to draw maps

The Data

• Need something to plot

Citibike Stations

• Station Data
• Map actual stations
• Count stations in each neighborhood

library(sf)
docks <- readr::read_csv(here::here('data', 'citibike_stations.csv')) |> 
    st_as_sf(
        # which columns to use as coordinates
        # always forget long/lat or lat/long
        coords=c('longitude', 'latitude'), 
        # keep the coordinate columns
        remove=FALSE,
        # projection system
        crs=4326
    )

• Pulled a few years ago
• Need to remember long/lat
• 4326 is lat/long

docks

Simple feature collection with 813 features and 18 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.09694 ymin: 40.6554 xmax: -73.90774 ymax: 40.81433
Geodetic CRS:  WGS 84
# A tibble: 813 × 19
      id stationName    availableDocks totalDocks latitude longitude statusValue
 * <dbl> <chr>                   <dbl>      <dbl>    <dbl>     <dbl> <chr>      
 1    72 W 52 St & 11 …             30         39     40.8     -74.0 In Service 
 2    79 Franklin St &…              1         33     40.7     -74.0 In Service 
 3    82 St James Pl &…             19         27     40.7     -74.0 In Service 
 4    83 Atlantic Ave …             35         62     40.7     -74.0 In Service 
 5   116 W 17 St & 8 A…              3         39     40.7     -74.0 In Service 
 6   119 Park Ave & St…             10         19     40.7     -74.0 In Service 
 7   120 Lexington Ave…              7         19     40.7     -74.0 In Service 
 8   127 Barrow St & H…             18         31     40.7     -74.0 In Service 
 9   128 MacDougal St …              0         30     40.7     -74.0 Not In Ser…
10   143 Clinton St & …              7         24     40.7     -74.0 In Service 
# … with 803 more rows, and 12 more variables: statusKey <dbl>,
#   availableBikes <dbl>, stAddress1 <chr>, stAddress2 <lgl>, city <lgl>,
#   postalCode <lgl>, location <lgl>, altitude <lgl>, testStation <lgl>,
#   lastCommunicationTime <chr>, landMark <lgl>, geometry <POINT [°]>

• sf object on top of tibble
• A lot of information about stations at a moment in time
• 813 features == rows
• 18 fields == columns
• Bounding box: rectangle that fits all data
• Geodetic CRS: lat/long

library(dplyr)
docks |> select(id, stationName, latitude, longitude)

Simple feature collection with 813 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.09694 ymin: 40.6554 xmax: -73.90774 ymax: 40.81433
Geodetic CRS:  WGS 84
# A tibble: 813 × 5
      id stationName                latitude longitude             geometry
   <dbl> <chr>                         <dbl>     <dbl>          <POINT [°]>
 1    72 W 52 St & 11 Ave               40.8     -74.0 (-73.99393 40.76727)
 2    79 Franklin St & W Broadway       40.7     -74.0 (-74.00667 40.71912)
 3    82 St James Pl & Pearl St         40.7     -74.0 (-74.00017 40.71117)
 4    83 Atlantic Ave & Fort Green…     40.7     -74.0 (-73.97632 40.68383)
 5   116 W 17 St & 8 Ave                40.7     -74.0  (-74.0015 40.74178)
 6   119 Park Ave & St Edwards St       40.7     -74.0 (-73.97803 40.69609)
 7   120 Lexington Ave & Classon A…     40.7     -74.0 (-73.95928 40.68677)
 8   127 Barrow St & Hudson St          40.7     -74.0 (-74.00674 40.73172)
 9   128 MacDougal St & Prince St       40.7     -74.0  (-74.00297 40.7271)
10   143 Clinton St & Joralemon St      40.7     -74.0  (-73.99338 40.6924)
# … with 803 more rows

• Ordinary {dplyr} functions
• The geometry column is sticky
  • Always selected even if not explicit
• Geometry matches long/lat

New York Neighborhoods

• Neighborhood boundaries
• Can be controversial

nyc <- read_sf(here::here('data', 'NYC-Neighborhood-Tabulation-Areas.geojson'))

• Geojson file
  • Like json, for geo
• Has geometry information
  • No need to build
• read_sf() can read many different types of files
  • database connections too

nyc |> select(-Shape__Area, -Shape__Length, -NTACode)

Simple feature collection with 195 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS:  WGS 84
# A tibble: 195 × 6
   OBJECTID BoroCode BoroName CountyFIPS NTAName                        geometry
      <int>    <int> <chr>    <chr>      <chr>                <MULTIPOLYGON [°]>
 1        1        3 Brooklyn 047        Borough Park  (((-73.97605 40.63128, -…
 2        2        4 Queens   081        Murray Hill   (((-73.80379 40.77562, -…
 3        3        4 Queens   081        East Elmhurst (((-73.8611 40.76367, -7…
 4        4        4 Queens   081        Hollis        (((-73.75726 40.71815, -…
 5        5        3 Brooklyn 047        Homecrest     (((-73.95859 40.61041, -…
 6        6        4 Queens   081        Fresh Meadow… (((-73.77758 40.7302, -7…
 7        7        4 Queens   081        St. Albans    (((-73.75205 40.70524, -…
 8        8        3 Brooklyn 047        Madison       (((-73.93754 40.60857, -…
 9        9        3 Brooklyn 047        Kensington-O… (((-73.97084 40.64639, -…
10       10        3 Brooklyn 047        Erasmus       (((-73.94827 40.63862, -…
# … with 185 more rows

• Multipolygon: Detached polygons make up one unit

Base/{sf}

• plot()

plot(docks)

• Plot for each column
• Color coded
• Named

docks |> select(geometry) |> plot()

• A little better
• Just the docks
• Nice to see position in city

nyc |> select(geometry) |> plot()
docks |> select(geometry) |> plot(col='blue', add=TRUE)

• Remember add=TRUE
• Some across the Hudson: Jersey City
• Can do more, but won’t

{ggplot2}

• Old standby
• Can plot anything

library(ggplot2)
ggplot(docks) + 
    geom_sf() + 
    coord_sf()

• Plots geometry column

• Station locations

• Can vaguely make out parts of NYC

• Nice to have the city plotted…

ggplot(docks) + 
    geom_sf(data=nyc) + 
    geom_sf() + 
    coord_sf()

• Sense for where docks are
• Looks like from the 80s
• Later on: better maps
• But first…

How Many Docks in Each Neighborhood?

• How many docs are in each neighborhood?

1. Find which docs are in which neighborhood
   • Spatial Join
2. Count the numbers in each neighborhood
   • Regular count
3. Plot the count
   • Choropleth

Spatial Join

• Find which docks in which neighborhoods

docks_in_neighborhoods <- docks |> 
    # only care about the station name and geometry
    select(stationName, geometry) |> 
    # spatial join
    st_join(
        # only need these columns from neighborhood tibble
        nyc |> select(BoroName, NTAName, geometry),
        # join rows where there is some over lap between
        # a dock and a neighborhod
        join=st_intersects,
        # keep only docks that are in a neighborhood
        # this will ignore Jersey City
        left=FALSE
    )

• st_join() joins rows where a dock (point) is in a neighborhood (polygon)
• Left join or inner join

docks_in_neighborhoods

Simple feature collection with 761 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.01713 ymin: 40.6554 xmax: -73.90774 ymax: 40.81433
Geodetic CRS:  WGS 84
# A tibble: 761 × 4
   stationName                               geometry BoroName  NTAName         
 * <chr>                                  <POINT [°]> <chr>     <chr>           
 1 W 52 St & 11 Ave              (-73.99393 40.76727) Manhattan Clinton         
 2 Franklin St & W Broadway      (-74.00667 40.71912) Manhattan SoHo-TriBeCa-Ci…
 3 St James Pl & Pearl St        (-74.00017 40.71117) Manhattan Chinatown       
 4 Atlantic Ave & Fort Greene Pl (-73.97632 40.68383) Brooklyn  Fort Greene     
 5 W 17 St & 8 Ave                (-74.0015 40.74178) Manhattan Hudson Yards-Ch…
 6 Park Ave & St Edwards St      (-73.97803 40.69609) Brooklyn  Fort Greene     
 7 Lexington Ave & Classon Ave   (-73.95928 40.68677) Brooklyn  Bedford         
 8 Barrow St & Hudson St         (-74.00674 40.73172) Manhattan West Village    
 9 MacDougal St & Prince St       (-74.00297 40.7271) Manhattan SoHo-TriBeCa-Ci…
10 Clinton St & Joralemon St      (-73.99338 40.6924) Brooklyn  Brooklyn Height…
# … with 751 more rows

• Shows which neighborhood each dock is in

count_by_neighborhood <- docks_in_neighborhoods |> 
    # remove geometry for fast counting
    st_drop_geometry() |> 
    count(NTAName) |> 
    # join the counts into the nyc neighborhood object
    right_join(nyc, by=c('NTAName'='NTAName')) |> 
    st_as_sf()

• Count number in each neighborhood
• Join back into neighborhood shapes
• Could have done differently
  • But this works for our purposes

count_by_neighborhood |> select(NTAName, n, BoroName, geometry)

Simple feature collection with 195 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS:  WGS 84
# A tibble: 195 × 4
   NTAName                                  n BoroName                  geometry
   <chr>                                <int> <chr>           <MULTIPOLYGON [°]>
 1 Astoria                                 22 Queens   (((-73.90203 40.76743, -…
 2 Battery Park City-Lower Manhattan       23 Manhatt… (((-74.00078 40.69429, -…
 3 Bedford                                 21 Brooklyn (((-73.94193 40.70073, -…
 4 Brooklyn Heights-Cobble Hill            11 Brooklyn (((-73.99237 40.6897, -7…
 5 Bushwick South                           3 Brooklyn (((-73.9323 40.70776, -7…
 6 Carroll Gardens-Columbia Street-Red…    34 Brooklyn (((-73.99237 40.6897, -7…
 7 Central Harlem North-Polo Grounds        4 Manhatt… (((-73.93445 40.83599, -…
 8 Central Harlem South                    12 Manhatt… (((-73.94177 40.8071, -7…
 9 Chinatown                               16 Manhatt… (((-73.98383 40.72148, -…
10 Clinton                                 14 Manhatt… (((-73.99383 40.77294, -…
# … with 185 more rows

• Have count along with neighborhood polygon

• One row per neighborhood

• Let’s plot it…

ggplot(count_by_neighborhood, aes(fill=n)) + 
    geom_sf() + 
    coord_sf()

• geom_sf() again
• Lighter blue: more docks
• Grey: no data
  • Jersey City is gone
• Other ways to map
• Both points and choropleth…

{tmap}

• Getting fancier

title

library(tmap)
tm_shape(docks) + 
    tm_dots()

• Pass data.frame to tm_shape()
• Draw with tm_dots()

tm_shape(nyc) +
    tm_polygons() +
    tm_shape(docks) + 
    tm_dots()

• Each data.frame goes to different tm_shape()
• Draw with tm_polygons() or tm_dots()

tm_shape(count_by_neighborhood) +
    tm_polygons(col='n')

• Column name as character

{leaflet}

• Interactive
• Plot docks…

library(leaflet)
leaflet(elementId='Docks') |>
    addTiles() |>
    addCircles(data=docks)

• elementId: for multiple plots in HTML
• addTiles(): Drawn right over Open Street Map
• Zoom and pan

leaflet(elementId='DocksHoods') |>
    addTiles() |>
    addPolygons(data=nyc, color='black', weight=3, opacity=0.2, popup=~NTAName) |>
    addCircles(data=docks, popup=~stationName)

• Different elementId
• Popup: shows name of dock or neighborhood on hover
• Tilde before column name
• Can choke on large data

library(leafgl)
leaflet(elementId='DocksHoodsGL') |>
    addTiles() |>
    addGlPolygons(data=nyc |> st_cast('POLYGON'), popup='NTAName', fillColor='black', fillOpacity=0.4, stroke=TRUE) |>
    addGlPoints(data=docks, popup='stationName')

• leafgl
• Uses WebGL
• addGlPolygons(): needs polygons not multipolygons: st_cast()
• Can’t see borders of polygons
• Column names in quotes
• Plot hundreds of thousands of points
• Choropleth…

dock_palette <- colorQuantile(palette="Dark2", domain=count_by_neighborhood$n, 
                              na.color="transparant")
leaflet(elementId='DocksHoodsCount') |> addTiles() |> 
    addPolygons(data=count_by_neighborhood, opacity=1, color='black',
                fillColor=~dock_palette(n), stroke=TRUE, weight=.5)

• Have to create the palette first
• Hard to see

leaflet(elementId='DocksHoodsCountZoomed') |> addTiles() |> 
    addPolygons(data=count_by_neighborhood, opacity=1, color='black',
                fillColor=~dock_palette(n), stroke=TRUE, weight=.5) |> 
    setView(lng=-73.85, lat=40.75, zoom=11)

• Default zoom: 10

• Zoom: 1-15

• A little better

• Try different base map

leaflet(elementId='DocksHoodsCountZoomedCarto') |> 
    addProviderTiles('CartoDB.Voyager', dependencyLocation=here::here()) |>
    addPolygons(data=count_by_neighborhood, opacity=1, color='black',
                fillColor=~dock_palette(n), stroke=TRUE, weight=.5) |> 
    setView(lng=-73.85, lat=40.75, zoom=11)

• Carto tiles
• One of many options
• Dev version: location of JavaScript files
• Not thirlled
• Try another way…

leaflet(elementId='DocksHoodsCountZoomedCartoGL') |> 
    addProviderTiles('CartoDB.Voyager', dependencyLocation=here::here()) |> 
    addGlPolygons(data=count_by_neighborhood |> st_cast("POLYGON"), 
                  opacity=1, fillColor='n') |> 
    setView(lng=-73.85, lat=40.75, zoom=11)

• Figures out colors for us
• Some areas that definitely don’t have bikes are colored
• Much easier to read
• Let’s add points…

leaflet(elementId='DocksHoodsCountZoomedCartoGLPoints') |> 
    addProviderTiles('CartoDB.Voyager', dependencyLocation=here::here()) |> 
    addGlPolygons(data=count_by_neighborhood |> st_cast("POLYGON"), 
                  opacity=1, fillColor='n') |> 
    addGlPoints(data=docks) |> 
    setView(lng=-73.85, lat=40.75, zoom=11)

• Docks definitely not in colored areas
• Looks like parks and wildlife refuges

Hex Bins

• Bin data into small geographic areas
• Density by area

library(h3jsr)
dock_hex <- docks |> 
    select(stationName) |> 
    mutate(hex=point_to_h3(geometry, res=8)) |> 
    st_drop_geometry() |> 
    count(hex) |> 
    h3_to_polygon(simple=FALSE)

dock_hex

Simple feature collection with 171 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -74.1013 ymin: 40.64706 xmax: -73.89869 ymax: 40.82119
Geodetic CRS:  WGS 84
First 10 features:
   n      h3_address h3_resolution                       geometry
1  2 882a100881fffff             8 POLYGON ((-73.97278 40.8022...
2  7 882a100883fffff             8 POLYGON ((-73.97033 40.7944...
3  1 882a100887fffff             8 POLYGON ((-73.98114 40.7967...
4  6 882a100889fffff             8 POLYGON ((-73.96443 40.8078...
5  8 882a10088bfffff             8 POLYGON ((-73.96198 40.7999...
6  6 882a100891fffff             8 POLYGON ((-73.95709 40.7841...
7  8 882a100893fffff             8 POLYGON ((-73.95464 40.7763...
8  3 882a100895fffff             8 POLYGON ((-73.96789 40.7865...
9  4 882a100897fffff             8 POLYGON ((-73.96544 40.7786...
10 9 882a100899fffff             8 POLYGON ((-73.94873 40.7897...

• H3: breaks world into hexagons
  • various sizes
• From Uber

leaflet(elementId='HexCounts') |> 
    addProviderTiles('CartoDB.Voyager', dependencyLocation=here::here()) |> 
    addGlPolygons(data=dock_hex, opacity=1, fillColor='n')

• Density plot by area
• Haven’t figured out legend with leafgl

Other Neat {leaflet} Features

• A lot of cool things we can do

library(leafem)
leaflet(elementId='CoolThings') |> 
    addProviderTiles('CartoDB.Voyager', dependencyLocation=here::here()) |> 
    addGlPolygons(data=dock_hex, opacity=1, fillColor='n') |> 
    addMeasure() |> 
    addMouseCoordinates() |> 
    addLogo(img=here::here('images', 'lander_analytics_sticker.png'), src='remote', url='https://www.landeranalytics.com')

• Measuring tools: upper right
• Mouse coordinates: upper left, can’t see here
• Logo: upper left

leaflet(elementId='MoreCoolThings') |> 
    addProviderTiles('CartoDB.Voyager', group='Streets', dependencyLocation=here::here()) |> 
    addProviderTiles('Esri.WorldImagery', group='Satellite', dependencyLocation=here::here()) |> 
    addGlPolygons(data=dock_hex, opacity=1, fillColor='n', group='Hexes') |>
    addCircles(data=docks, group='Docks') |> addMeasure() |> addMouseCoordinates() |> 
    addLogo(img=here::here('images', 'lander_analytics_sticker.png'), src='remote', url='https://www.landeranalytics.com') |> 
    addLayersControl(overlayGroups=c('Docks', 'Hexes'), baseGroups=c("Satellite", "Open Streets"),
                     options=layersControlOptions(collapsed=FALSE))

• Switch between maps
• Show hex, points or both

Game Time

• Use GIS to solve a game
• Most popular game…

• Not this game
• But …

• Map version of Wordle
• Excuse the poor resolution

• Map version of Wordle
• Outline of country
• After guess
  • How far from correct country
  • Direction to correct country
• No clue what this is

guesses <- tibble::tribble(
    ~Country, ~Distance, ~Direction,
    'Iceland', 13427, 'South',
    'Sierra Leone', 7144, 'South',
    'Lesotho', 3404, 'Southwest'
)

guesses

# A tibble: 3 × 3
  Country      Distance Direction
  <chr>           <dbl> <chr>    
1 Iceland         13427 South    
2 Sierra Leone     7144 South    
3 Lesotho          3404 Southwest

• Distance: Center of guess to center of correct country
• Had a feeling a small island somewhere
• Let’s find the centers of those countries

data(countries110, package='rnaturalearth')

starting <- countries110 |> 
    # this is an sp object so we make it into sf
    st_as_sf() |> 
    # make sure it's lat/long
    st_transform(4326) |> 
    select(brk_name) |> 
    # this is an sp object so we make it into sf
    inner_join(guesses, by=c('brk_name'='Country'))

# this finds the centroids of each country
centers <- starting |> st_centroid()

centers

Simple feature collection with 3 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -18.76554 ymin: -29.62479 xmax: 28.17182 ymax: 65.07986
Geodetic CRS:  WGS 84
      brk_name Distance Direction                   geometry
1      Iceland    13427     South POINT (-18.76554 65.07986)
2      Lesotho     3404 Southwest POINT (28.17182 -29.62479)
3 Sierra Leone     7144     South POINT (-11.79541 8.529459)

• World map from Natural Earth
• Find centers: st_centroid()

leaflet(elementId='WorldleGuesses') |> 
    addTiles() |> 
    addPolygons(data=starting, fillColor='green', stroke=FALSE, opacity=1) |> 
    addMarkers(data=centers)

• Markers at centers of countries
• Polygons don’t line up perfectly

circles <- centers |> 
    # we use the distance from each center
    # this is stored in km so we multiply by 1000 to get meters
    st_buffer(dist=centers$Distance*1000) |> 
    # get just the outline of the cirles
    st_cast("LINESTRING")

circles

Simple feature collection with 3 features and 3 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -180 ymin: -60.55612 xmax: 178.7879 ymax: 73.37476
Geodetic CRS:  WGS 84
      brk_name Distance Direction                       geometry
1      Iceland    13427     South LINESTRING (-36.46156 -54.9...
2      Lesotho     3404 Southwest LINESTRING (-5.850577 -45.4...
3 Sierra Leone     7144     South LINESTRING (-45 -48.94435, ...

• Compute circle of radius distance from each center
• Using spherical geometry: S2 from Google
• LINESTRING: Gives the outline of the circle

leaflet(elementId='WorldleCircles') |> 
    addTiles() |> 
    # we use the color column defined earlier
    addGlPolylines(data=circles, popup='brk_name', color='brk_name') |> 
    addMarkers(data=centers, label=~brk_name)

• Big distances
• Wonky near poles
• Purple: weird shape because Iceland’s circle goes on other side of globe
• Triangulate where three circle intercept
  • Flat purple part doesn’t count

overlaps <-circles |> 
    st_intersection() |> 
    filter(n.overlaps >= 3)

overlaps

Simple feature collection with 1 feature and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3.483787 ymin: -54.73521 xmax: 3.483787 ymax: -54.73521
Geodetic CRS:  WGS 84
    brk_name Distance Direction n.overlaps origins                   geometry
1.2  Iceland    13427     South          3 1, 2, 3 POINT (3.483787 -54.73521)

• Find where they intersect: st_intersection()
• Triangulation: 3 circles overlap

leaflet(elementId='Triangulation') |> 
    addTiles() |>
    addCircles(data=overlaps) |> 
    # 100 km search area
    addPolylines(data=overlaps |> st_buffer(dist=100*1000))

• Bouvette Island
• So small needed to search around triangulation
• Not in Natural Earth

leaflet(elementId='Winner') |> 
    addProviderTiles('Esri.WorldImagery', dependencyLocation=here::here()) |> 
    setView(lng=3.3558, lat=-54.42852, zoom=12)

• Uninhabited nature reserve

Conclusions

• Wrapping Up

• sf objects are tibbles with geometry information
• {sf} package offers most standard GIS operations
  • Others in {lwgeom}, {h3jsr}, {geodist}
• Use st_join() for spatial joins
• Make maps with
  • {ggplot2}
  • {tmap}
  • {leaflet}
• All about points, lines and polygons

• sf objects are tibbles with geometry information
• sf package offers most standard GIS operations
• Use st_join() for spatial joins
• Make maps with
  • ggplot2
  • tmap
  • leaflet
• All about points, lines and polygons

Thank You

• Thank You