Skip to content

Instantly share code, notes, and snippets.

@adamhsparks
Forked from andrewheiss/geom_point_sf.md
Created August 20, 2018 20:16
Show Gist options
  • Save adamhsparks/f44f26c0fb630a03d1f1401502a74882 to your computer and use it in GitHub Desktop.
Save adamhsparks/f44f26c0fb630a03d1f1401502a74882 to your computer and use it in GitHub Desktop.
geom_point and geom_sf
library(tidyverse)
library(sf)

Generate fake data (tribble() is a cool function, btw)

events <- tribble(
  ~event, ~latitude, ~longitude,
  "event 1", 41.46, -74.53,
  "event 2", 43.59, -89.79,
  "event 3", 40.49, -80.11,
  "event 4", 32.77, -83.65
)

Latitude and longitude are just numbers

events

#> # A tibble: 4 x 3
#>     event latitude longitude
#>     <chr>    <dbl>     <dbl>
#> 1 event 1    41.46    -74.53
#> 2 event 2    43.59    -89.79
#> 3 event 3    40.49    -80.11
#> 4 event 4    32.77    -83.65

Use st_as_sf() to convert latitude and longitude into the magic geometry column. Use 4326 as the CRS, since it's the World Geodetic System, the standard -180 to 180 lat/long system used by GPS systems

events_sf <- events %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

events_sf

#> Simple feature collection with 4 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -89.79 ymin: 32.77 xmax: -74.53 ymax: 43.59
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#> # A tibble: 4 x 2
#>     event          geometry
#>     <chr>  <simple_feature>
#> 1 event 1 <POINT (-74.5...>
#> 2 event 2 <POINT (-89.7...>
#> 3 event 3 <POINT (-80.1...>
#> 4 event 4 <POINT (-83.6...>

Now it'll plot with geom_sf()

ggplot() + geom_sf(data = events_sf)

BUUUUT. There are issues with geom_sf() and points. Like, if you try to resize them with size = 3, barely anything happens. So you'll want to use geom_point() instead. You're still using other geom_sf stuff, though, like the shapefiles, and it's helpful to convert the lat/lon numbers into an sf geometry object, because you can transform the coordinates to other projections, like the Albers projection for the contiguous US (102003)

So transform the geometry, etc., and then extract the points with st_coordinates() and join those points back to the main data frame, then use the extracted points in geom_point()

events_transformed <- events_sf %>%
  st_transform(crs = 102003)

# The coordinates are different now!
events_transformed

#> Simple feature collection with 4 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 498908.3 ymin: -454520 xmax: 1764805 ymax: 698100.5
#> epsg (SRID):    102003
#> proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
#> # A tibble: 4 x 2
#>     event          geometry
#>     <chr>  <simple_feature>
#> 1 event 1 <POINT (17648...>
#> 2 event 2 <POINT (49890...>
#> 3 event 3 <POINT (13292...>
#> 4 event 4 <POINT (11470...>

# st_coordinates() extracts the lon/lat values as a data frame with X and Y columns
st_coordinates(events_transformed)

#>           X         Y
#> 1 1764804.7  643907.8
#> 2  498908.3  698100.5
#> 3 1329260.8  446475.9
#> 4 1147033.0 -454520.0

# Bind the coordinates to the full data frame
events_transformed_with_lat_lon <- cbind(events_transformed, st_coordinates(events_transformed))

# The X and Y points are now extracted!
events_transformed_with_lat_lon

#> Simple feature collection with 4 features and 3 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 498908.3 ymin: -454520 xmax: 1764805 ymax: 698100.5
#> epsg (SRID):    102003
#> proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
#>     event         X         Y                       geometry
#> 1 event 1 1764804.7  643907.8 POINT (1764804.73588095 643...
#> 2 event 2  498908.3  698100.5 POINT (498908.343714951 698...
#> 3 event 3 1329260.8  446475.9 POINT (1329260.75728861 446...
#> 4 event 4 1147033.0 -454520.0 POINT (1147032.98687954 -45...

It works!

ggplot() +
  # geom_sf(data = whatever_shapefile_youre_using) +
  geom_point(data = events_transformed_with_lat_lon, aes(x = X, y = Y), size = 4) +
  coord_sf(crs = 102003)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment