Nubian Square Median Household Income

library(tigris)
library(acs)
library(leaflet)
library(mapview)
library(ggplot2)
library(maptools)
library(stringr)

shapefile <- counties(state='MA', class='sp')

## Warning in proj4string(obj): CRS object has comment, which is lost in output

plot(shapefile)

head(shapefile)

## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : -73.06851, -69.85886, 41.18705, 42.72161  (xmin, xmax, ymin, ymax)

## Warning in proj4string(x): CRS object has comment, which is lost in output

## crs         : +proj=longlat +datum=NAD83 +no_defs 
## variables   : 17
## names       : STATEFP, COUNTYFP, COUNTYNS, GEOID,       NAME,          NAMELSAD, LSAD, CLASSFP, MTFCC, CSAFP, CBSAFP, METDIVFP, FUNCSTAT,      ALAND,     AWATER, ... 
## min values  :      25,      001, 00606927, 25001, Barnstable, Barnstable County,   06,      H1, G4020,   148,  12700,    14454,        A, 1021019039, 1004283854, ... 
## max values  :      25,      027, 00606940, 25027,  Worcester,  Worcester County,   06,      H4, G4020,   148,  49340,    14454,        N, 3912554985,  666897066, ...

plot(shapefile)
pointLabel(coordinates(shapefile),labels=shapefile$COUNTYFP)

countylist <- c('17','21','25', '09')
shapefile <- tracts(state='MA', county=countylist, year=2018, class='sp')

## Warning in proj4string(obj): CRS object has comment, which is lost in output

plot(shapefile)

# Span
myspan <- 5

# End year
myendyear <- 2018

# Which table
mytable <- acs.lookup(endyear=2018, table.number="B19013")

## Warning in acs.lookup(endyear = 2018, table.number = "B19013"): acs.lookup for endyear>2015: using 2015 variable codes to access 2018 data.
##   (See ?acs.lookup for details)

## Warning in acs.lookup(endyear = 2018, table.number = "B19013"): temporarily downloading and using archived XML variable lookup files;
##   since this is *much* slower, recommend running
##   acs.tables.install()

# Limit the geography
# Convert the list of counties that is in text form to numeric form
countylist2 <- as.numeric(countylist)
# Define the geography to limit the data import
mygeo <- geo.make(state=25, county=countylist2, tract="*")

results(mytable)$variable.name

## [1] "Median household income in the past 12 months (in 2015 Inflation-adjusted dollars)"

myvars <- mytable[1]

mydata <- acs.fetch(endyear=myendyear, span=myspan, geography=mygeo, variable=myvars)

mystate <- as.character(mydata@geography$state)

mycounty <- as.character(mydata@geography$county)

# add leading zeros to make the total number of characters 3
mycountypadded <- str_pad(mycounty, 3, pad="0")

mytract <- as.character(mydata@geography$tract)

acsgeoid <- paste0(mystate,mycountypadded,mytract)

# Put the dataset into a dataframe format
mydatadf <- data.frame(acsgeoid, mydata@estimate)



# Rename column names
colnames(mydatadf)=c("GEOID", "medianincome")


# There is some weirdness with some tracts showing up as 100%, so eliminate those by subsetting the results to exclude them.
mydatadf <- subset(mydatadf, mydatadf$medianincome>0)

# Merge the ACS data with the shapefile
mydatamerged <- geo_join(shapefile, mydatadf, "GEOID", "GEOID")

mypopup <- paste0("GEOID: ", mydatamerged$GEOID, "<br>", "Median Household Income: ", round(mydatamerged$medianincome,1), "%")

mypal <- colorNumeric(
  palette = "YlGnBu",
  domain = NULL
)

# Store values into variables
myLAT <- 42.366740
myLNG <- -71.084760
myZOOM <- 12
myTILES <- "CartoDB.Positron"

mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = myZOOM) %>%
  addPolygons(data = mydatamerged, 
              fillColor = ~mypal(medianincome), # this makes the fill color on the spectrum of the palette
              color = "white", # this is the color of the outline of the shapes
              fillOpacity = 0.7, # how see through the shapes are
              weight = 1, 
              smoothFactor = 0.2,
              popup = mypopup) %>%
  addLegend(pal = mypal, 
            values = mydatamerged$medianincome, 
            position = "bottomright", 
            title = "Median Household Income",
            labFormat = labelFormat(suffix = "%"),
            bins=5)

# Draw the map
mymap