This lab complements the exercise using Qgis for today

Here, we use tidycensus to read some tract data, learn its projection information, transform it to a new coordinate system and measure some distance between features.

library(tidycensus)
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Read in Bexar county tracts

sa_acs<-get_acs(geography = "tract",
                state="TX",
                county = c("Bexar"),
                year = 2017,
                variables=c( "DP05_0001E", 
                            "DP03_0119PE") ,
                geometry = T, output = "wide")
## Getting data from the 2013-2017 5-year ACS
## Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## Using the ACS Data Profile
#create a county FIPS code - 5 digit
sa_acs$county<-substr(sa_acs$GEOID, 1, 5)

#rename variables and filter missing cases
sa_acs2<-sa_acs%>%
  mutate(totpop= DP05_0001E, ppov=DP03_0119PE) %>%
#  st_transform(crs = 102740)%>%
  na.omit()

find coordinate system of current map

st_crs(sa_acs2)
## Coordinate Reference System:
##   EPSG: 4269 
##   proj4string: "+proj=longlat +datum=NAD83 +no_defs"

create basic map

ppov_map<-sa_acs2 %>%
  mutate(cpov=cut(ppov,breaks = quantile(ppov, na.rm=T, p=seq(0,1,length.out = 6)),include.lowest = T))

library(ggsn)
## Loading required package: ggplot2
## Loading required package: grid
p1<-ggplot(ppov_map, aes(fill = cpov, color = cpov)) + 
  geom_sf() + 
  ggtitle("Proportion in poverty", 
          subtitle = "Bexar County Texas, 2017 - Quantile Breaks")+
    scale_fill_brewer(palette = "Reds") + 
  scale_color_brewer(palette = "Reds")+
  theme(axis.text.x = element_blank(), axis.text.y = element_blank())+
  north(ppov_map)+
  scalebar(ppov_map, location="bottomleft", dist=5, transform = T,dist_unit = "km",  model="WGS84", st.size =2 )
p1

re-project map into South Central Texas projection

new_sa<-st_transform(sa_acs2, crs = 2278)

#Extract two tracts
twtr<-new_sa%>%
  filter(GEOID %in% c(48029181820, 48029110600))

# get centroid coordinates for two tracts
tr_co<-st_centroid(twtr)
## Warning in st_centroid.sf(twtr): st_centroid assumes attributes are constant
## over geometries of x
#Measure feet apart
st_distance(tr_co)
## Units: [US_survey_foot]
##          [,1]     [,2]
## [1,]     0.00 64043.26
## [2,] 64043.26     0.00