POST
|
After some messing I solved this with minimal use of arcpy. I used arcpy in the first code block to create the pickled dictionary referenced in the previous post and then the second code block I use to very quickly calculate higher order neighbors. By removing arcpy from the mix computation time goes down from hours to seconds. Designed to work with a TIGER street file. Hopefully, this will save someone the the days of suffering I just endured! To create a dictionary of first order street network neighbors: ##This Script Creates and pickles a dictionary of street connectivity
##Designed to work with a TIGER street file
########################
##IMPORT LIBRARIES
########################
import cPickle as pickle
import arcpy
#Set workspace to geodatabase
arcpy.env.workspace = "aGDB.gdb/aFC"
arcpy.env.overwriteOutput = True
##################
##CREATE VARIABLES
##################
st2neighs = {}
#Make layers
arcpy.MakeFeatureLayer_management("SSOStreetsMerged", "strLyr")
#Make cursor of all streets
sts = arcpy.SearchCursor("SSOStreetsMerged", "", "", "", "")
#Create a dictionary containing each street and its neighbors
for ast in sts:
st = int(ast.getValue("TLID"))
if st not in st2neighs:
st2neighs[st] = []
arcpy.SelectLayerByAttribute_management("strLyr", "NEW_SELECTION", "\"TLID\" = " + str(st))
arcpy.SelectLayerByLocation_management("strLyr", "INTERSECT", "strLyr", "", "NEW_SELECTION")
tsts = arcpy.SearchCursor("strLyr", "", "", "", "")
for tst in tsts:
if int(tst.getValue("TLID")) != st:
if int(tst.getValue("TLID")) not in st2neighs[st]:
st2neighs[st].append(int(tst.getValue("TLID")))
##Pickle seg2st dictionary
pfile = open("st2neighs.pkl", "wb")
pickle.dump(st2neighs, pfile)
pfile.close() Then to trace the higher order network neighbors without using arcpy (select by location): ##Select K-order neighbors from a set of sampled streets.
##saves in dictionary format such that
##the key is the sampled street and the neighboring streets are the values
##################
##IMPORT LIBRARIES
##################
import random as random
import cPickle as pickle
import copy
#######################
##LOAD PICKLED DATA
#######################
seg_file = open("seg2st.pkl", "rb")
st_file = open("st2neighs.pkl", "rb")
seg2st = pickle.load(seg_file)
st2neighs = pickle.load(st_file)
st_file.close()
seg_file.close()
st2neigh = copy.deepcopy(st2neighs)
##################
##FUNCTIIONS
##################
##Takes in a dict of segments (key) and their streets (values).
##returns the desired number of sampled streets per segment
##returns a dict keyed segment containing tlids.
def selectSample(seg2st, nbirths):
randSt = {}
for segK in seg2st.iterkeys():
ranSamp = [int(random.choice(seg2st[segK])) for i in xrange(nbirths)]
randSt[segK] = []
for aSamp in ranSamp:
z = copy.deepcopy(aSamp)
randSt[segK].append(z)
return randSt
##Takes in a list of all streets (keys) and their first order neighbors (values)
##Takes in a list of sampled streets
##returns a dict of all sampled streets and their neighbors.
##Higher order selections should be possible with findMoreNeighbors
##logic is the same but replacing sample (input) with output from
##findFirstNeighbors
def initList(sample):
compSts = {}
for samp in sample:
for rSt in sample[samp]:
if rSt not in compSts:
compSts[rSt] = []
return compSts
def findFirstNeighbors(st2neigh, compSts):
for x in compSts:
st = copy.deepcopy(st2neigh )
for z in st:
compSts .append(z)
return compSts
def findMoreNeighbors(st2neigh, compSts):
for aSt in compSts:
stx = copy.deepcopy(compSts[aSt])
for st in stx:
nStx = copy.deepcopy(st2neigh[st])
for nSt in nStx:
if nSt not in compSts[aSt]:
compSts[aSt].append(nSt)
return compSts
#####################
##ANALYSIS
#####################
samp = selectSample(seg2st, 1)
n0 = initList(samp)
n1 = findFirstNeighbors(st2neigh, copy.deepcopy(n0))
n2 = findMoreNeighbors(st2neigh, copy.deepcopy(n1))
n3 = findMoreNeighbors(st2neigh, copy.deepcopy(n2))
n4 = findMoreNeighbors(st2neigh, copy.deepcopy(n3))
n5 = findMoreNeighbors(st2neigh, copy.deepcopy(n4))
n6 = findMoreNeighbors(st2neigh, copy.deepcopy(n5))
n7 = findMoreNeighbors(st2neigh, copy.deepcopy(n6))
n8 = findMoreNeighbors(st2neigh, copy.deepcopy(n7))
n9 = findMoreNeighbors(st2neigh, copy.deepcopy(n8))
n10 = findMoreNeighbors(st2neigh, copy.deepcopy(n9))
#####################
##CHECK RESULTS
#####################
def checkResults(neighList):
cntr = {}
for c in neighList.iterkeys():
cntr = 0
for a in neighList :
cntr += 1
return cntr
c1 = checkResults(n1)
c2 = checkResults(n2)
c3 = checkResults(n3)
c4 = checkResults(n4)
c5 = checkResults(n5)
c6 = checkResults(n6)
c7 = checkResults(n7)
... View more
07-27-2011
09:01 PM
|
1
|
0
|
182
|
POST
|
I would like to identify the Kth order neighbors of an edge on a network, specifically the neighbors of a large set of streets. For example, I have a street that I'm interested in looking at, call this the focal street. For each focal street I want to find the streets that share an intersection, these are the first order neighbors. Then for each of those streets that share an intersection with the focal street I would like to find their neighbors (these would be the second order neighbors), and so on... Calculating the first order neighbors using ArcGIS' geoprocessing library (arcpy) took 6+ hours, second order neighbors has been running for 20+ hours. Needless to say I want to find a more efficient solution. I have created a python dictionary which is keyed on each street and contains the connected streets as values. For example: st2neighs = {street1: [street2, street3, street5], street2: [street1, street4], ...}. street 1 is connected to street 2, 3, 5; street 2 is connected to street 1 and 4; etc.. There are around 12000 streets in the study area most have fewer than 7 connected streets. I have tried to do this using python-only and have not succeeded (see here). SelectLayerByLocation_management seems insanely slow. My read of the forums has turned up little, one person suggests recursively buffering each street and then making an in-memory layer but it is unclear to me how this would add efficiency (it increases the required number SelectLayerByLocation calls, doesn't it?). A pickled version of the st2neighs used in the code below is here. The intended behavior is: For each street in the city (focal street): [INDENT]For each of the focal streets neighbors (focal street neighbor): [INDENT]Select the neighbors of the focal street neighbor using SelectLayerByLocation [/INDENT][/INDENT] Any thoughts on why the code below is taking so long to run? Did I make a mistake? By my estimates the code involves ~ 72000 calls to SelectLayerByLocation. The code has been running about 20 hours so that's about 1 second per SelectLayerByLocation (and it isn't done yet). Is this normal? Can the efficiency be improved somehow? Thanks for your help and code follows... ##This Script Creates and pickles a dictionary of street connectivity
########################
##IMPORT LIBRARIES
########################
import cPickle as pickle
import arcpy
##################
##LOAD DATA
##################
#Set workspace to geodatabase
arcpy.env.workspace = "aGDB.gdb/aLayer"
arcpy.env.overwriteOutput = True
st_file = open("st2neighs.pkl", "rb")
#st2neighs contains first order neighbors
st2neighs = pickle.load(st_file)
#Make layers
arcpy.MakeFeatureLayer_management("SSOStreetsMerged", "strLyr")
#Second order neighbors
for street in st2neighs:
for astreet in st2neighs[street]:
arcpy.SelectLayerByAttribute_management("strLyr", "NEW_SELECTION", "\"ID\" = " + str(astreet))
arcpy.SelectLayerByLocation_management("strLyr", "INTERSECT", "strLyr", "", "NEW_SELECTION")
neighs = arcpy.SearchCursor("strLyr", "", "", "", "")
for aneigh in neighs:
neigh = int(aneigh.getValue("ID"))
if neigh != street:
if neigh not in st2neighs[street]:
st2neighs[street].append(neigh)
pfile = open("st2neighs2.pkl", "wb")
pickle.dump(st2neighs, pfile)
pfile.close()
... View more
07-27-2011
05:41 PM
|
0
|
2
|
1792
|
Title | Kudos | Posted |
---|---|---|
1 | 07-27-2011 09:01 PM |
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|