# shapefile sample --------------------------------------------------------
library(rgdal)
library(sf)
library(ggplot2)
library(tidyverse)
library(tigris)
library(leaflet)
# data source -------------------------------------------------------------
# <https://www.census.gov/geo/maps-data/data/cbf/cbf_counties.html>
# <https://www.census.gov/geo/maps-data/data/cbf/cbf_state.html>
# use rgdal ---------------------------------------------------------------
# <https://www.nceas.ucsb.edu/scicomp/usecases/ReadWriteESRIShapeFiles>
# optionally report shapefile details
ogrInfo("./input/shape/census", "cb_2016_us_county_500k")
ogrInfo("./input/shape/census", "cb_2016_us_state_500k")
# read in shapefiles
counties <- readOGR("./input/shape/census", "cb_2016_us_county_500k")
states <- readOGR("./input/shape/census", "cb_2016_us_state_500k")
# note that readOGR will read the .prj file if it exists
print(proj4string(counties))
print(proj4string(states))
# generate a simple map showing all three layers
plot(counties, axes=TRUE, border="gray")
plot(states, axes=TRUE, border="gray")
# write out a new shapefile (including .prj component)
writeOGR(counties, "./input/shape/census", "counties-rgdal", driver="ESRI Shapefile")
# use sf package ----------------------------------------------------------
counties_sf <- st_read("./input/shape/census", "cb_2016_us_county_500k")
states_sf <- st_read("./input/shape/census", "cb_2016_us_state_500k")
# add county centroids
counties_sf <- counties_sf %>%
mutate(
centroid_x = as.numeric(st_coordinates(st_geometry(st_centroid(.)))[,"X"]),
centroid_y = as.numeric(st_coordinates(st_geometry(st_centroid(.)))[,"Y"])
)
# add state names to county df, filter to Colorado
county_names <- counties_sf %>%
left_join(as.data.frame(states_sf),
by = "STATEFP",
suffix = c(".cty", ".state")
) %>%
select(STATEFP,
NAME.cty,
NAME.state,
STUSPS,
centroid_x,
centroid_y,
geometry.cty) %>%
rename(geometry = geometry.cty)
# plot with ggplot
p <- county_names %>% filter(NAME.state == "Colorado") %>% ggplot +
theme_bw() +
geom_sf(color = "black", fill = "light grey") +
geom_text(
check_overlap = TRUE,
size = 3,
aes(x = .$centroid_x,
y = .$centroid_y,
label = .$NAME.cty)
) +
labs(
title = paste0("State of Colorado")
)
# use tigris package to download shape files ------------------------------
library(tigris)
library(leaflet)
tigris_states <- states(cb = TRUE, class = "sf")
leaflet(states) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons(fillColor = "white",
color = "black",
weight = 0.5) %>%
setView(-98.5795, 39.8282, zoom=3)
# when you have shp, dbf, prj, and shx ------------------------------------
library(sf)
my_sf <- st_read("shapetest.shp")
ggplot() + geom_sf(data = my_sf, alpha = 0.2, aes(fill = MUKEY_V)) +
coord_sf(datum = st_crs(my_sf)$epsg)