R/nyc_gis_adjacency.R

In davewutchiett/nets_demo_spatial_dynamics: Predicting Change and Modeling Relationships In NYC Local Retail and Population.

#' get adjacency matrix from GIS files to define network structure
library(maptools)
library(spdep)
wd = ""
df.file = ""
setwd(wd)
df = read.csv(df.file)
id = 'trtid10'
alt_id = 'GEOID10'

no_nb_tracts = c(36005051600, 36081091601, 36081107201,36005000100) #exclude island tracts without neighbors
df = df[! df[,id] %in% no_nb_tracts,]
valid_ids = (unique(df[,id]))
ny_map_file = "ny_gis/tl_2010_36_tract10.shp"
nj_map_file = "nj_gis/tl_2010_34_tract10.shp"
ct_map_file = "ct_gis/tl_2010_09_tract10.shp"
pa_map_file = "pa_gis/tl_2010_42_tract10.shp"
ny_map <- readShapePoly(ny_map_file, IDvar = alt_id)
nj_map <- readShapePoly(nj_map_file, IDvar = alt_id) 
ct_map <- readShapePoly(ct_map_file, IDvar = alt_id)
pa_map <- readShapePoly(pa_map_file, IDvar = alt_id)
maps = spRbind( spRbind(spRbind( ny_map, nj_map), ct_map), pa_map)
maps = maps[maps$GEOID10 %in% valid_ids, ]
# plot(maps) # view map of NYC metro Census tracts
nb.lists <- poly2nb(maps, row.names = maps$GEOID10 , queen = F) #neighbors if share common edges with length >0
adj_mat <- as.data.frame(nb2mat(nb.lists, style = "B", zero.policy = TRUE)) #to 'n x n' adjacency matrix with binary indicators
colnames(adj_mat) = rownames(adj_mat)
write.csv(adj_mat, 'nyc_adjacency_matrix.csv')