alathrop
12/20/2017 - 8:28 PM

shape files


# 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)