How close is everything?

This file provides example code used for the following study of parcel-level accessibility to businesses in Long Beach, Califonria following this paper in the Professional Geographer and this Metropolitan Futures Initiative report .

This script uses sample data for blocks and business establishments in Irvine, California. The distance from each block to each of 10 (fake) types of business establishments (denoted by the field bus_type) is measured in two ways: the distance to the nearest establishment of that type, and a measure of all establishments within one mile.

The network analyst output - which is a line shapefile for each establishment type - is then joined back to the original blocks shapefile in R. The nearest establishment measure generates two attributes: distance (in feet) and drive time (in minutes) to the nearest establishment. The count-within-a-mile measure generates three attributes: the count of establishments within 1 mile, a gravity measure using inverse distance decay, and a gravity measure using inverse squared distance decay. See the report and paper above for details.

This tutorial contains four sections. Run the preamble and file declaration codeblocks, then one of these:

(1) Run in ArcGIS. Find the nearest establishment of each type.

(2) Run in R. Nearest establishment: join Network Analyst output to block shapefile.

(3) Run in ArcGIS. Find all establishments within 1 mile.

(4) Run in R. Join count of establishments in 1 mile and gravity measure to block shapefile.



In [1]:

    
# Preamble to use for all the below scripts. Make sure to set your working directory as 'wkspace':
import arcpy
import numpy as np
import datetime
from arcpy.da import*
from arcpy.sa import*
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Network")
def unique_values(table, field):
    with arcpy.da.SearchCursor(table, [field]) as cursor:
        return sorted({row[0] for row in cursor})
wkspace = "C:\\Users\\Research\\Project"
arcpy.env.workspace = wkspace









    Out[1]:





u'CheckedOut'



In [ ]:

    
#### (1) Find The Nearest Establishment of each type #### 
print 'start: ', datetime.datetime.now()

# Read in files.  Note: network dataset is NOT found in GitHub Repository.
nd = "\\your_network_dataset"
establishments = "\\Irvine_businesses_randomized.shp"
arcpy.MakeFeatureLayer_management(establishments, "estabs_lyr")
types = unique_values(establishments, "bus_type")
blocks = "\\Irvine_blocks.shp"
arcpy.MakeFeatureLayer_management(blocks, "blocks_lyr")

# Run network analyst for each establishment type 
for i in types:
    arcpy.SelectLayerByAttribute_management("estabs_lyr", "NEW_SELECTION", " \"bus_type\" = '%s' "%i)
    outNALayer = arcpy.na.MakeODCostMatrixLayer(nd, "od_lyr", "Distance", "#", 1, ["Distance", "Minutes"])
    arcpy.na.AddLocations("od_lyr", "Origins", "blocks_lyr", "#", "5000 Meters")
    arcpy.na.AddLocations("od_lyr", "Destinations", "estabs_lyr", "#", "5000 Meters")
    arcpy.na.Solve("od_lyr", "SKIP", "CONTINUE")
    subLayerNames = arcpy.na.GetNAClassNames(outNALayer.getOutput(0))
    lines = subLayerNames["ODLines"]
    arcpy.FeatureClassToFeatureClass_conversion(lines, wkspace, "\\lines_near_%s.txt"%i)

print 'end: ', datetime.datetime.now()



In [ ]:

    
#### (2) Nearest establishment: join Network Analyst output to block shapefile in R

setwd("C:/Users/Research/Project")
library(foreign)

bl = read.dbf("Irvine_blocks.dbf")
bl$ESRI_na = as.numeric(bl$FIDnum + 1)
estabs = read.dbf("Irvine_businesses_randomized.dbf")
types = as.character(unique(unlist(estabs$bus_type))) 
types = types[!is.na(types)]
types = types[!types %in% c("none")]

for(i in 1:length(types)){
  out = read.dbf(paste("./lines_near_", types[i], ".dbf", sep=""))
  bl[,ncol(bl)+1] = out$Total_Dist[match(bl$ESRI_na, out$OriginID)]
  colnames(bl)[ncol(bl)] = paste("dist_", types[i], sep="")
  bl[,ncol(bl)+1] = out$Total_Minu[match(bl$ESRI_na, out$OriginID)]
  colnames(bl)[ncol(bl)] = paste("mins_", types[i], sep="")
}

write.dbf(bl, "Irvine_blocks.dbf")



In [4]:

    
#### (3) Find all establishments within 1 mile. 
print 'start: ', datetime.datetime.now()

# Read in files.  Note: network dataset is NOT found in GitHub Repository.
nd = "\\your_network_dataset"
establishments = "\\Irvine_businesses_randomized.shp"
arcpy.MakeFeatureLayer_management(establishments, "estabs_lyr")
types = unique_values(establishments, "bus_type")
blocks = "\\Irvine_blocks.shp"
arcpy.MakeFeatureLayer_management(blocks, "blocks_lyr")

# Run network analyst for each establishment type 
for i in types:
    arcpy.SelectLayerByAttribute_management("estabs_lyr", "NEW_SELECTION", " \"bus_type\" = '%s' "%i)
    outNALayer = arcpy.na.MakeODCostMatrixLayer(nd, "od_lyr", "Distance", 5280, "", ["Distance", "Minutes"])
    arcpy.na.AddLocations("od_lyr", "Origins", "blocks_lyr", "#", "5000 Meters")
    arcpy.na.AddLocations("od_lyr", "Destinations", "estabs_lyr", "#", "5000 Meters")
    arcpy.na.Solve("od_lyr", "SKIP", "CONTINUE")
    subLayerNames = arcpy.na.GetNAClassNames(outNALayer.getOutput(0))
    lines = subLayerNames["ODLines"]
    arcpy.FeatureClassToFeatureClass_conversion(lines, wkspace, "lines_1mi_%s"%i) 

print 'finish ', datetime.datetime.now()



In [ ]:

    
#### (4) Join count of establishments in 1 mile and gravity measure to block shapefile using R. 

setwd("C:/Users/Research/Project")
library(foreign)
library(doBy)

bl = read.dbf("Irvine_blocks.dbf")
bl$ESRI_na = as.numeric(bl$FIDnum + 1)
estabs = read.dbf("Irvine_businesses_randomized.dbf")
types = as.character(unique(unlist(estabs$bus_type))) 
types = types[!is.na(types)]
types = types[!types %in% c("none")]

for(i in types){
  out = read.dbf(paste(".lines_1mi_", types[i], ".dbf", sep=""))
  out$invd_mile = 1/((out$Total_Dist+0.02)/5280)
  out$invsq_mile = 1/((out$Total_Dist+0.02)/5280)^2
  out$est = 1
  outeach = summaryBy(invd_mile + invsq_mile + est ~ OriginID, FUN=sum, data=out)
  bl[,ncol(parcels)+1] = outeach$invd_mile[match(bl$ESRI_na, outeach$OriginID)]
  bl[,ncol(bl)][is.na(bl[,ncol(bl)])] = 0
  colnames(bl)[ncol(bl)] = paste("invmi_", types[i], sep="")
  bl[,ncol(bl)+1] = outeach$invsq_mile[match(bl$ESRI_na, outeach$OriginID)]
  bl[,ncol(bl)][is.na(bl[,ncol(bl)])] = 0
  colnames(bl)[ncol(bl)] = paste("invsq_", types[i], sep="")
  bl[,ncol(bl)+1] = outeach$est.sum[match(bl$ESRI_na, outeach$OriginID)]
  bl[,ncol(bl)][is.na(bl[,ncol(bl)])] = 0
  colnames(bl)[ncol(bl)] = paste("num_est_", types[i], sep="")
}

write.dbf(bl, "Irvine_blocks.dbf")



In [ ]: