Day: November 11, 2023💾

You can do IDW interpolation of missing Census Tracts fairly easily in R using the gstat library. The key is to make sure you use a projected dataset. Other interpolation methods are covered here: https://rspatial.org/raster/analysis/4-interpolation.html

library(tidycensus)
library(tidyverse)
library(raster)
library(gstat)

# obtain census data on veteran status by tract and then
# reproject the shapefile geometery into a projected coordinate system

acs <- get_acs("tract",
               state='ny',
               survey='acs5', 
               var=    
                 c('Total'='B21001_001',
                   'Veteran'='B21001_002'
                 ),
               cache_table = T,
               geometry = T, 
               resolution='20m',
               year = 2020, 
               output = "wide"
) %>% st_transform(26918)

# calculate the percentage of veterans per census tract
acs <- mutate(acs, vet_per = VeteranE/TotalE)

# create a copy of census tracts, dropping any NA values
# from vet_per field
vetNA <- acs %>% drop_na(vet_per)

# a raster should be created to do interpolation into
r <- raster(vetNA, res=1000)


# set the foruma based on field (vet_per) that contains
# the veterans percent to interpolate. This use IDW interpolation
# for all points, weighting farther ones less
gs <- gstat(formula = vet_per~1, locations = vetNA)

# interpolate the data (average based on nearby points)
nn <- interpolate(r, gs)

# extract the median value of the raster interpolation from the original shapefile,
# into a new column set as est
acs<- cbind(acs, est=exactextractr::exact_extract(nn, acs, fun='median'))

# replace any NA values with interpolated data so the map doesn't contain
# holes. You should probably mention that missing data was interpolated when
# sharing your map.
acs <- acs %>% mutate(vet_per = ifelse(is.na(vet_per), est, vet_per))