Working with Maps and GIS data
Notable packages
The tidycensus packages
- Website: https://walkerke.github.io/tidycensus/
- Install with
install.packages("tidycensus") - Vignettes: https://walkerke.github.io/tidycensus/articles/basic-usage.html
- Reference: https://walkerke.github.io/tidycensus/reference/index.html
The simple feature packages
- Website: https://github.com/r-spatial/sf
- Install with
install.packages("sf")- On Mac, you may need to install
udunitsfirst, for example, with homebrew - On Linux, install
libudunits2-devwith your software management system first
- On Mac, you may need to install
- Vignettes: https://r-spatial.github.io/sf/articles/sf1.html
- Reference: https://r-spatial.github.io/sf/reference/index.html
Code examples
Data files used in the examples
- Portland Metro 1994 TAZ shape file: Download & unzip to
datasubdirectory of your RStudio project folder - Portland Metro 1994 Survey geocode.raw.zip
Map ACS 2012-16 Median Household Income by Census Tract
## Install the tidycensus package if you haven't yet
#install.packages("tidycensus")
library(tidycensus)
library(ggplot2)
library(dplyr)
## setup cenus api key
## signup your census api key at http://api.census.gov/data/key_signup.html
#census_api_key("YOUR_CENSUS_API_KEY") #
portland_tract_medhhinc <- get_acs(geography = "tract",
year = 2016, # 2012-2016
variables = "B19013_001", # Median Household Income in the Past 12 Months
state = "OR",
county = c("Multnomah County", "Washington County", "Clackamas County"),
geometry = TRUE) # load geometry/gis info## Getting data from the 2012-2016 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 12%
|
|======== | 13%
|
|========= | 15%
|
|========== | 16%
|
|=========== | 17%
|
|============ | 18%
|
|============= | 20%
|
|============== | 21%
|
|=============== | 23%
|
|================ | 24%
|
|================= | 26%
|
|================= | 27%
|
|================== | 28%
|
|=================== | 29%
|
|==================== | 31%
|
|===================== | 32%
|
|====================== | 34%
|
|======================= | 35%
|
|======================== | 36%
|
|======================== | 38%
|
|========================= | 39%
|
|========================== | 40%
|
|=========================== | 42%
|
|============================ | 43%
|
|============================= | 45%
|
|============================== | 46%
|
|=============================== | 47%
|
|================================ | 49%
|
|================================= | 50%
|
|================================= | 51%
|
|================================== | 53%
|
|=================================== | 54%
|
|==================================== | 56%
|
|===================================== | 57%
|
|====================================== | 58%
|
|======================================= | 59%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 64%
|
|========================================== | 65%
|
|=========================================== | 66%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================== | 70%
|
|=============================================== | 72%
|
|================================================ | 73%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 77%
|
|=================================================== | 79%
|
|==================================================== | 80%
|
|===================================================== | 81%
|
|====================================================== | 83%
|
|====================================================== | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|========================================================= | 88%
|
|========================================================== | 90%
|
|=========================================================== | 91%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================== | 95%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|=================================================================| 99%
|
|=================================================================| 100%(myplot <- ggplot(portland_tract_medhhinc) +
geom_sf(aes(fill = estimate)) +
coord_sf(datum = NA) + theme_minimal())
ggsave("output/mymap.pdf", myplot)## Saving 7 x 5 in imageInteractive Maps of ACS 2012-16 Median Household Income by Census Tract
## Install the mapview package if you haven't yet
#install.packages("mapview")
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.2library(mapview)
library(dplyr)
mapview(portland_tract_medhhinc %>% select(estimate),
col.regions = sf.colors(10), alpha = 0.1)Example of spatial analysis: spatial join
library(sf)
library(readr)
# read 1994 Metro TAZ shape file
taz_sf <- st_read("data/taz1260.shp", crs=2913)## Reading layer `taz1260' from data source `/home/lmwang/datascience2018/data/taz1260.shp' using driver `ESRI Shapefile'## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that## Simple feature collection with 1247 features and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 7435706 ymin: 447796.5 xmax: 7904748 ymax: 877394.8
## epsg (SRID): 2913
## proj4string: +proj=lcc +lat_1=46 +lat_2=44.33333333333334 +lat_0=43.66666666666666 +lon_0=-120.5 +x_0=2500000.0001424 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=ft +no_defs# read geocode.raw file that contains X and Y coordinates
portland94_df <- read_csv("data/portland94_geocode.raw.zip", col_names=c("uid", "X", "Y", "case_id",
"freq", "rtz", "sid",
"totemp94", "retemp94"))## Multiple files in zip: reading 'geocode.raw'## Parsed with column specification:
## cols(
## uid = col_double(),
## X = col_double(),
## Y = col_double(),
## case_id = col_integer(),
## freq = col_integer(),
## rtz = col_integer(),
## sid = col_integer(),
## totemp94 = col_double(),
## retemp94 = col_double()
## )portland94_df <- portland94_df %>%
filter(X!=0, Y!=0) %>%
sample_n(500)
# create a point geometry with x and y coordinates in the data frame
portland94_sf <- st_as_sf(portland94_df, coords = c("X", "Y"), crs = 2913)
# spatial join to get TAZ for observations in portland94_sf
portland94_sf <- st_join(portland94_sf, taz_sf)
head(portland94_sf)## Simple feature collection with 6 features and 8 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 7588418 ymin: 655251.6 xmax: 7677074 ymax: 731410
## epsg (SRID): 2913
## proj4string: +proj=lcc +lat_1=46 +lat_2=44.33333333333334 +lat_0=43.66666666666666 +lon_0=-120.5 +x_0=2500000.0001424 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=ft +no_defs
## uid case_id freq rtz sid totemp94 retemp94 TAZ
## 1 2063744104 1475 48 216 1 1596.065 270.1489 216
## 2 2030982205 14557 1 940 1 7716.001 1059.1665 940
## 3 5057563203 12284 43 1053 2 8277.523 313.1616 1053
## 4 2143082103 17159 1 890 1 4150.559 1092.1593 890
## 5 2073613204 19969 32 606 1 5022.285 772.7363 606
## 6 2085892107 11959 25 793 1 16642.590 2620.7116 793
## geometry
## 1 POINT (7588418 680452.1)
## 2 POINT (7655822 655251.6)
## 3 POINT (7650362 731410)
## 4 POINT (7663721 691092)
## 5 POINT (7677074 694609)
## 6 POINT (7649698 675806.8)ggplot() +
geom_sf(data = taz_sf, aes(alpha=0.9)) +
geom_sf(data = portland94_sf)