Calculating state mean centers in R

meancenterweights

library(tidyverse)
library(sf)

#Join population data to tract boundaries and convert to point coordinates
tractdata<-read_csv("nhgis0105_2017_tract_racepov.csv") #Data on population downloaded from NHGIS
tracts<-st_read("nhgis0101_shapefile_tl2017_us_tract_2017/US_tract_2017.shp",stringsAsFactors=FALSE) %>% #Tract boundaries downloaded from NHGIS
  st_transform(4326) 
tract_center<-tracts %>%
  st_centroid() %>%
  st_coordinates() %>%
  as_tibble() %>%
  cbind(tracts %>% st_set_geometry(NULL) %>% select(GISJOIN)) 

#Calculate weighted mean centers by states
tract_mean_center<-tractdata %>%
  filter(AHZAE001>0) %>% #Total population variable
  left_join(tract_center) %>% 
  group_by(STATE) %>%
  mutate(Xw=X*AHZAE001,
         Yw=Y*AHZAE001) %>%
  filter(is.na(Xw)==FALSE) %>%
  summarise(meanX=sum(Xw)/sum(AHZAE001),
            meanY=sum(Yw)/sum(AHZAE001)) %>%
  st_as_sf(coords=c("meanX","meanY"),crs=4326,remove=FALSE)

#Load place names and calculate closest one to each state center
place_points<-st_read("nhgis0110_shapefile_tl2017_us_place_2017/US_place_2017.shp", #Place boundaries downloaded from NHGIS
                      stringsAsFactors=FALSE) %>%
  st_transform(4326) %>%
  st_centroid() %>%
  rownames_to_column(var = "value") %>%
  mutate(value=as.numeric(value))
matches<-st_nearest_feature(tract_mean_center,place_points) %>%
  as_tibble() %>%
  bind_cols(tract_mean_center %>% st_set_geometry(NULL)) %>%
  left_join(place_points %>% 
              st_set_geometry(NULL) %>%
              select(NAME,GEOID,value)) %>%
  st_as_sf(coords=c("meanX","meanY"),crs=4326)

st_write(matches,"state_centers.gpkg")