I'd add my answer as a comment, but since I can't comment (no reputation enough), I'll add as an answer.
It seems that you are trying to add centroid connectors to feed into a traffic assignment algorithm. Right? I should say that doing it automatically is not considered best practice in this industry, but it you are going to do it, I'd like to suggest you to make it proportional to the amount of network comprised within the polygon (roadway stock) and the displacement of this network as well.
I'll also mix some of the transportation modelling jargon with the GIS jargon...
First of all, are you considering splitting links or just connecting the new centroid connectors (dummy links) to the existing extremities of the links?
Since I don't know what is the case, I'll assume that you want to connect them to the existing extremities. And it DOES make a huge difference in the algorithm if anybody is wondering.
Now is where I think there are two major questions you need to think about:
1st - How to weigh the criteria closest links versus busiest links
2nd - How to spread the connectors through the TAZ so you don't have them all clustered together?
I'll not answer thse questions, so when it comes to the part of the algorithm to actually add centroids, I'll add a single one.
The algorithm I am suggesting is:
- Compute the number of connectors to be added in each polygon (for each TAZ)
- Associate this number of connectors with the centroids file
- Extract all nodes from the road network (link extremities) to a new layer
- Tag TAZ (polygon) ID to each one of these nodes
Loop through the centroids adding the connectors to a new layer
Now let's put some codes to these steps (I hope the variables names are self explanatory):
Steps 1 & 2:
Let's assume there is TAZ field in both the centroids layer and the TAZ (polygon) layer.
Something like this:
connec={}
features = taz_layer.getFeatures()
idx=taz_layer.fieldNameIndex("TAZ")
for f in features:
taz_id=feat.attributes()[idx]
geom = f.geometry()
area= geom.area()
connectors=area*number_of_connectors_per_unit_of_area
connec[taz_id]=connectors
features = centroids_layer.getFeatures()
idx=centroids_layer.fieldNameIndex("TAZ")
for f in features:
taz_id=feat.attributes()[idx]
a=centroids_layer.dataProvider().changeAttributeValues({f.id(): {idx:connec[taz_id]}})
centroids_layer.updateFields()
Step 3:
#Add node layer as suggested in the PyQGIS cookbook
writer = QgsVectorFileWriter("node_layer.shp", "CP1250", fields, \
QGis.WKBPoint, None, "ESRI Shapefile")
idx=taz_layer.fieldNameIndex("TAZ")
for feat in line_layer.getFeatures():
link=list(feat.geometry().asPolyline())
node_a=(link[0][0],link[0][1])
node_b=(link[-1][0],link[-1][1])
new_feat = QgsFeature()
new_feat.setGeometry(QgsGeometry.fromPoint(QgsPoint(node_a)))
new_feat.setAttributes(["how busy is the link", "name of the street", "etc", "etc])"#Here you need to capture and add all the characteristics of the links you want to consider when adding connectors.
writer.addFeature(new_feat)
Step 4:
#This code is very slow, but will suffice
#You should look into using spatial index, but this is what I could write by heart...
idx=taz_layer.fieldNameIndex("TAZ")
idx2=node_layer.fieldNameIndex("TAZ")
for feat in node_layer.getFeatures(): #The extremities
for feature in taz_layer.getFeatures():
enclosed = feature.geometry().contains(feat.geom())
if enclosed:
taz_id=feature.attributes()[idx]
a=node_layer.dataProvider().changeAttributeValues({feat.id(): {idx2:taz_id}})
centroids_layer.updateFields()
Step 5:
#Add link layer as suggested in the PyQGIS cookbook
writer = QgsVectorFileWriter("node_layer.shp", "CP1250", fields, \
QGis.WKBLinestring, None, "ESRI Shapefile")
index = QgsSpatialIndex()
for feature in node_layer.values():
index.insertFeature(feature)
for feature in centroids_layer.values():
geom=pt.feature.geometry()
pt = QgsPoint(geom.asPoint())
nearest = index_from.nearestNeighbor(pt, 100) #100 is a guess, you might need to do more...
#Nearest contains the IDS of the nodes, so it is wise to build dictionaries or numpy arrays with all the relevant criteria (TAZ, if it is busy or not, etc) before you even start looping through centroids
for i in nearest:
if criteria1(i) and criteria2(i) and ... and criterian(i):
P=[pt, node_i.geometry().asPoint()]
line_geom=QgsGeometry.fromPolyline(P)
new_feat = QgsFeature()
new_feat.setGeometry(line_geom)
new_feat.setAttributes['Whatever','attributes','you','want']
writer.addFeatures([new_feat])
Well... It ended up being an answer longer than I had anticipated, but it should get you started...
And just as a warning, I did NOT test this code, so it is likely to have typos and other mistakes...
Best Answer
dassouki, can you work with Hawth's Analysis Tools for ArcGIS 9.x?
---EDIT---
with "Add XY to Table (points)" (in "Table tools" from Hawth's bar), create the fields with the coordinates from each point from the point layer.
in excel, for example, open the csv file and add 4 fields: fromX, fromY, toX, toY and copy the coordinates from each point (work with filters, is quite fast) and then save the csv file
in arcmap, use "Add XY line data from table (in "Table tools" from Hawth's bar) in a line layer using the csv file as source and the various connections are created
in the line layer simbology, use "Multiple attributes - Quantity by category", using Weight field
this is the result of my test: