library(tidycensus)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(dplyr)
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
sa_acs$county <- substr(sa_acs$GEOID, 1, 5)
sa_acs2 <- sa_acs %>%
mutate(
totpop = DP05_0001E,
ppov = DP03_0119PE
) %>%
st_transform(crs = 4269) %>%
na.omit()
st_crs(sa_acs2)
#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: grid
p1 <- ggplot(ppov_map, aes(fill= cpov, color = cpov)) +
geom_sf() +
ggtitle("Proportion in poverty",
subtitle = "Bexar county Texas, 2017 - Quantile breaks/ NAD83") +
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
## Distance using NAD83
#new_sa <- st_transform(sa_acs2, crs = 2278)
twtr <- sa_acs2 %>%
filter(GEOID %in% c(48029181820, 48029110600))
st_distance((twtr))
## Units: [m]
## [,1] [,2]
## [1,] 0.00 16646.58
## [2,] 16646.58 0.00
sa_acs3 <- sa_acs %>%
mutate(
totpop = DP05_0001E,
ppov = DP03_0119PE
) %>%
st_transform(crs = 3083) %>%
na.omit()
st_crs(sa_acs3)
#create basic map
ppov_map3 <- sa_acs3 %>%
mutate(
cpov = cut(ppov, breaks = quantile(ppov, na.rm = T, p = seq(0,1, length.out = 6)), include.lowest = T)
)
p2 <- ggplot(ppov_map3, aes(fill= cpov, color = cpov)) +
geom_sf() +
ggtitle("Proportion in poverty",
subtitle = "Bexar county Texas, 2017 - Quantile breaks / Texa Centric") +
scale_fill_brewer(palette = "Reds") +
scale_color_brewer(palette = "Reds") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank()) +
north(ppov_map)
p2
## Distance using NAD83 / Texas Centric Albers Equal Area
#new_sa <- st_transform(sa_acs2, crs = 2278)
twtr <- sa_acs3 %>%
filter(GEOID %in% c(48029181820, 48029110600))
st_distance((twtr))
## Units: [m]
## [,1] [,2]
## [1,] 0.0 16661.4
## [2,] 16661.4 0.0
Looks like they’re different by
x <- 16661.4 - 16646.58
print(x)
## [1] 14.82
meters
In general, you need, not only an accurate system of projection, but also a consistant system of projection; there would be errors between layers for the same point if they were on two different projections. The method used to measure distance should also be in a useful unit like meters or miles, not something hard to comprehend like radians.