Last active
April 6, 2025 11:55
-
-
Save zross/857978baba866f7894cc3fbe2fa52868 to your computer and use it in GitHub Desktop.
Spatial analysis with R
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# An example of spatial analysis with R that involves grabbing data | |
# from a PostGIS database, joining to tabular data, projecting | |
# simplifying and creating an interactive map. | |
# This data is publicly available but requires some processing | |
# if you want the county boundaries you can find them here: | |
# cftp://ftp2.census.gov/geo/tiger/TIGER2016/COUNTY/ | |
# The tabular data requires some processing but the raw data | |
# is here: | |
# https://www.census.gov/data/datasets/2016/demo/saipe/2016-state-and-county.html | |
# If you want to learn more about processing data with the sf | |
# package you can go to this DataCamp course: | |
# https://www.datacamp.com/courses/spatial-analysis-in-r-with-sf-and-raster | |
library(dplyr) | |
library(sf) | |
library(RColorBrewer) | |
library(mapview) | |
# Read the poverty data | |
pov <- readr::read_csv("data.csv") | |
# Make a connection to a local PostGIS database | |
con<-RPostgreSQL::dbConnect(RPostgreSQL::PostgreSQL(), | |
host="localhost", | |
user= "postgres", | |
dbname="data_census") | |
# Read data from PostGIS to R with sf pacakge | |
cnty <- st_read_db(con, table = c("geo", "tl_2016_us_county")) | |
cnty <- mutate(cnty, geoid = paste0(statefp, countyfp)) | |
cnty <- st_transform(cnty, crs = 2163) # project | |
# Original size of file | |
object.size(cnty) # 127 megabytes | |
# Simplify | |
cnty_simp <- st_simplify(cnty,dTolerance = 500, preserveTopology = TRUE) | |
# New size of file | |
object.size(cnty_simp) # 5 megabytes | |
# Link to tabular data | |
cnty_dat <- filter(pov, year == 2016) %>% | |
inner_join(cnty_simp, ., by = "geoid") | |
cnty_dat <- filter(cnty_dat, !substring(geoid,1,2)%in%c("02", "15")) | |
# Create an interactive map | |
mapviewOptions(legend.pos = "bottomright") | |
mapview::mapview(cnty_dat, zcol="pov_per", | |
col.regions = colorRampPalette(brewer.pal(7, "YlOrBr")), | |
legend = TRUE, | |
layer.name = "Percent poverty", | |
lwd = 0.2) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment