library(tidycensus)
library(tidyverse)
library(sf)
library(ggplot2)
library(classInt)
kz_acs <- get_acs(geography = "tract",
                state="MI",
                county = "Kalamazoo",
                year = 2019,
                variables=c("DP05_0001E", "DP03_0119PE", "DP03_0062E") ,
                geometry = T,
                output = "wide")
## Getting data from the 2015-2019 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
head(kz_acs)
# create a county FIPS code - 5 digit
kz_acs$county<-substr(kz_acs$GEOID, 1, 5)

kz_acs2<-kz_acs%>%
  mutate(totpop= DP05_0001E,
         ppov=DP03_0119PE,
         hhmed=DP03_0062E) %>%
  st_transform(crs = 2919)%>%
  na.omit()

kz_acs$county<-substr(kz_acs$GEOID, 1, 5)

kz_acs2<-kz_acs%>%
  mutate(totpop= DP05_0001E,
         ppov=DP03_0119PE,
         hhmed=DP03_0062E) %>%
  st_transform(crs = 2253)%>%
  na.omit()
# create a county FIPS code - 5 digit
kz_acs$county<-substr(kz_acs$GEOID, 1, 5)

kz_acs3<-kz_acs%>%
  mutate(totpop= DP05_0001E,
         ppov=DP03_0119PE,
         hhmed=DP03_0062E) %>%
  st_transform(crs = 3032)%>%
  na.omit()

kz_acs4<-kz_acs%>%
  mutate(totpop= DP05_0001E,
         ppov=DP03_0119PE,
         hhmed=DP03_0062E) %>%
  st_transform(crs = 2278)%>%
  na.omit()
st_crs(kz_acs2)
## Coordinate Reference System:
##   User input: EPSG:2253 
##   wkt:
## PROJCRS["NAD83 / Michigan South (ft)",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Michigan South zone (International feet)",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",41.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-84.3666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",43.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",42.1,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",13123359.58,
##             LENGTHUNIT["foot",0.3048],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["foot",0.3048],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["foot",0.3048]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["foot",0.3048]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United States (USA) - Michigan - counties of Allegan; Barry; Bay; Berrien; Branch; Calhoun; Cass; Clinton; Eaton; Genesee; Gratiot; Hillsdale; Huron; Ingham; Ionia; Isabella; Jackson; Kalamazoo; Kent; Lapeer; Lenawee; Livingston; Macomb; Mecosta; Midland; Monroe; Montcalm; Muskegon; Newaygo; Oakland; Oceana; Ottawa; Saginaw; Sanilac; Shiawassee; St Clair; St Joseph; Tuscola; Van Buren; Washtenaw; Wayne."],
##         BBOX[41.69,-87.2,44.22,-82.13]],
##     ID["EPSG",2253]]
library(tmap)
library(tmaptools)

tm_shape(kz_acs2)+
  tm_polygons("ppov", title="Household Median Income", palette="Blues", style="quantile", n=5 )+
  tm_format("World", title="Kalamazoo Household Median Income - Quantile Breaks \n Using South Michigan CRS", legend.outside=T)+
  tm_scale_bar()+
  tm_compass()

library(tmap)
library(tmaptools)

tm_shape(kz_acs3)+
  tm_polygons("ppov", title="Household Median Income", palette="Blues", style="quantile", n=5 )+
  tm_format("World", title="Kalamazoo Household Median Income - Quantile Breaks \n Using Australia CRS", legend.outside=T)+
  tm_scale_bar()+
  tm_compass()

library(tmap)
library(tmaptools)

tm_shape(kz_acs4)+
  tm_polygons("ppov", title="Household Median Income", palette="Blues", style="quantile", n=5 )+
  tm_format("World", title="Kalamazoo Household Median Income - Quantile Breaks \n Using South Texas CRS", legend.outside=T)+
  tm_scale_bar()+
  tm_compass()

new_kz <- st_transform(kz_acs2, crs=2253)

#Extract two tracts
twtr<-new_kz%>%
  filter(GEOID %in% c(26077002101, 26077002004))

# 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: [ft]
##          [,1]     [,2]
## [1,]     0.00 18450.73
## [2,] 18450.73     0.00
aus_kz <- st_transform(kz_acs2, crs=3032)

austwtr<-aus_kz%>%
  filter(GEOID %in% c(26077002101, 26077002004))

# get centroid coordinates for two tracts
aus_tr_co<-st_centroid(austwtr)
## Warning in st_centroid.sf(austwtr): st_centroid assumes attributes are constant
## over geometries of x
#Measure feet apart
st_distance(aus_tr_co)
## Units: [m]
##          [,1]     [,2]
## [1,]     0.00 33014.27
## [2,] 33014.27     0.00
tx_kz <- st_transform(kz_acs2, crs=2278)

txtwtr<-tx_kz%>%
  filter(GEOID %in% c(26077002101, 26077002004))

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

The distances between the two centroids are different depending on which coordinate system is used.