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)
Since this was posted in 2017, the {sf} package has become more particular with how CRS systems are defined. You now need to use
st_crs()
and provide a string of the CRS system, prefaced with the abbreviation for the catalog (ESRI, EPSG, etc.). See here for more details.Here's how the gist above works with
st_crs()
: