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.7.2, GDAL 2.4.2, PROJ 5.2.0
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) %>%
  na.omit()

Find Coordinate System of Current Map

st_crs(sa_acs2)
## Coordinate Reference System:
##   EPSG: 4269 
##   proj4string: "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

Create Basic Map

library(classInt)
ppov_map<-sa_acs2 %>%
  mutate(cpov=cut(ppov,breaks = quantile(ppov, na.rm=T, p=seq(0,1,length.out = 6)),include.lowest = T),
jpov = cut(ppov,breaks=data.frame(classIntervals(var=sa_acs2$ppov, n=5, style="pretty")[2])[,1], include.lowest = T))
library(ggsn)
## Loading required package: ggplot2
## Loading required package: grid
library(RColorBrewer)

p1<-ggplot(ppov_map, aes(fill = cpov)) +
  #geom_sf()+
  geom_sf(aes(fill = cpov), color = "black", size=0.05) + 
  ggtitle("Proportion of Population in Poverty", 
          subtitle = "Bexar County Texas, 2017 - Quantile Breaks")+
    scale_fill_brewer(name="Percent in Poverty", palette = "Reds") + 
    #scale_color_brewer(name="", palette = "Greys") + 
  #theme_void()+
  theme(plot.title=element_text(size = 14, face = "bold"), axis.text.x = element_blank(), axis.title.x = element_blank(), axis.text.y = element_blank(), axis.title.y = element_blank(), axis.ticks=element_blank(), plot.background=element_blank(), panel.background = element_blank()) +
  north(ppov_map) +
  scalebar(ppov_map, location="bottomleft", dist=5, transform = T,dist_unit = "km",  model="WGS84", st.size =2)

p2<-ggplot(ppov_map, aes(fill = jpov)) +
  #geom_sf()+
  geom_sf(aes(fill = jpov), color = "black", size=0.1) + 
  ggtitle("Proportion of Population in Poverty", 
          subtitle = "Bexar County Texas, 2017 - Pretty Breaks") +
    scale_fill_brewer(name="Percent in Poverty", palette = "PuBu") + 
    #scale_color_brewer(name="", palette = "Greys") + 
  #theme_void() +
  theme(plot.title=element_text(size = 14, face = "bold"), axis.text.x = element_blank(), axis.title.x = element_blank(), axis.text.y = element_blank(), axis.title.y = element_blank(), axis.ticks=element_blank(), plot.background=element_blank(), panel.background = element_blank()) +
  north(ppov_map)+
  scalebar(ppov_map, location="bottomleft", dist=5, transform = T,dist_unit = "km",  model="WGS84", st.size =2)
p1

p2

Re-Project Map Into Texas South Central Projection

new_sa<-st_transform(sa_acs, crs = 2278)
plot(new_sa)

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

# 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
plot(tr_co)

#Measure feet apart
st_distance(tr_co)
## Units: [US_survey_foot]
##          [,1]     [,2]
## [1,]     0.00 64043.26
## [2,] 64043.26     0.00
SouthCentDistMiles<-st_distance(tr_co)/1609.344 #miles
SouthCentDistMiles
## Units: [US_survey_foot]
##          [,1]     [,2]
## [1,]  0.00000 39.79463
## [2,] 39.79463  0.00000

Re-Project Map Into Texas Centric Albers Equal Area Projection

new_sa2<-st_transform(sa_acs2, crs = 3083)
plot(new_sa2)

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

plot(twtr2)

# get centroid coordinates for two tracts
tr_co2<-st_centroid(twtr2)
## Warning in st_centroid.sf(twtr2): st_centroid assumes attributes are
## constant over geometries of x
plot(tr_co2)

#Measure meters apart
st_distance(tr_co2)
## Units: [m]
##          [,1]     [,2]
## [1,]     0.00 19536.81
## [2,] 19536.81     0.00
AlbersDistMiles<-st_distance(tr_co2)/1609.344 #miles
AlbersDistMiles
## Units: [m]
##          [,1]     [,2]
## [1,]  0.00000 12.13961
## [2,] 12.13961  0.00000

Response

What is the distance between the two points?

39.79 miles using the Texas South Center coordinate reference system (crs=2278)

Is this distance interpretable?

The distance is in meters (64043.26) and unit conversion is necessary to see distance in miles (39.79).

How does it compare to the one you got using the Texas South Central projection?

The final reproject is using the Texas Centric Albers Equal Area coordinate reference system (crs=3083) and produced different results. This distance was less at 12.14 miles compared to 39.79 miles with the Texas South Central coordinate reference system.

In general, why is it important to have an accurate system of projection? How could your results be sensitive to this?

It is important to have an accurate system of projection to produce appropriate, interpretable, and useable results. Results are sensitive and dependent on the coordinate reference system or geographical/spatial frame of reference. Output and measures depend on a consistent frame of reference across the entire geographical area of interest. Disparate coordinate reference system with geographic area will produce distorted results and with greater distance between points of interest, distortions will be more pronounced.