# -*- coding: utf-8 -*- """ Created on Thu Jan 7 12:07:26 2021 Geopreprocessing script This script takes inputs (neighbourhood household points,road network, neighbourhood exit point, transfer station and depot points) and uses ArcGIS python modules to create a graph of nodes and links. It also outputs required CSV files of the characteristics of the graph for later processing. It is written to be run cell-by-cell. @author: Charlie Ferguson """ #%% Import Packages import arcpy as ap import numpy as np import pandas as pd import os import progressbar #%% Input Data Location="Lima" #Neighbourhood Identifier SpaCRS="WGS 1984 UTM Zone 18S" #WGS1984 UTM Co-ordinate Reference System HouseholdCentroidPoints="{HouseholdsShapefile}.shp" #Point shapefile [Attribute Table: FID Shape Id] RoadNetwork_Full="{RoadNetworkShapefile}.shp"#Polyline shapefile ExitPoint="{ExitPointShapefile}.shp" #Point shapefile [Attribute Table: FID Shape Id] TransferStations="{TransferStationsShapefile}.shp" #Point shapefile [Attribute Table: FID Shape Id] RoadVerticesTransferStations="{RoadNetworkVerticesandTransferStationLocationsShapefile}.shp" #Point shapefile [Attribute Table: FID Shape Id] RoadVerticesTransferStations_RV_Num=199 #Number of vertices in the road network RoadVerticesTransferStations_TS_Num=222 #Number of vertices for transfer stations #### All input files should be in a file entitled 'GeoData' within the root directory (see below) #%% Define primary variales NumUsers=10 #Number of CBS users in the neighbourhood NumSets=50 #Number of random sets of CBS users #%% Define file connections os.chdir("{Path to root folder}") FilePath=os.path.abspath(os.getcwd()) #Define File path GeometricOutput_Folder=FilePath+"\\GeometricOutput\\"+Location+"\\" #Output file for geometric data #%% Processing ap.env.workspace=FilePath+"\\GeoData\\"+Location #Define Arcpy workspace ap.AddXY_management(HouseholdCentroidPoints) #Ensure xy data in shapefile #To run this for loop the folders 'Contributing shapefiles" should be empty bar = progressbar.ProgressBar(maxval=NumSets).start() for i in range(0,NumSets): #i=0 ap.Delete_management(["RoadVerticesTransferStations_EraseBuffer.shp","RoadNetwork_TransferVerticesErased.shp", "RoadProximityCBSUsers_PointsTable.csv","RoadProximityCBSUsers_PointsShape.shp", "RoadProximityCBSUsers_PointsShape_ConflictFree.shp","ExitPointA_Buffer.shp", "RoadVerticesTransferStations_withoutExit.shp","Final_Points.shp","Final_Lines.shp", "FinalProximityTable.csv" ]) #### #Create Random CBS Users (and save) ap.management.CreateRandomPoints(GeometricOutput_Folder, "CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i), constraining_feature_class=HouseholdCentroidPoints, number_of_points_or_field=NumUsers) #### #Erase all stationary points from road network ap.Buffer_analysis(RoadVerticesTransferStations, "RoadVerticesTransferStations_EraseBuffer.shp","1 Meter") ap.Erase_analysis(in_features=RoadNetwork_Full, erase_features="RoadVerticesTransferStations_EraseBuffer.shp", out_feature_class="RoadNetwork_TransferVerticesErased.shp") RoadNetwork='RoadNetwork_TransferVerticesErased.shp' ## #### #Identify Closest Road Points ap.GenerateNearTable_analysis(GeometricOutput_Folder+"\\CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i)+".shp",RoadNetwork,"RoadProximityCBSUsers_PointsTable.csv",location="LOCATION",closest = 'CLOSEST') ap.XYTableToPoint_management(in_table="RoadProximityCBSUsers_PointsTable.csv", out_feature_class="RoadProximityCBSUsers_PointsShape.shp", x_field="NEAR_X", y_field="NEAR_Y",coordinate_system=ap.SpatialReference(SpaCRS)) DF_walkingdistances = pd.read_csv(FilePath+"\\GeoData\\"+Location+"\\RoadProximityCBSUsers_PointsTable.csv") DF_walkingdistances=DF_walkingdistances[["IN_FID","NEAR_DIST"]] np.savetxt(FilePath+"\\ConstructedCSVFiles\\"+Location+"\\"+ "CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i) + "_WalkingDistances.csv",DF_walkingdistances, delimiter=",") del DF_walkingdistances ap.Integrate_management("RoadProximityCBSUsers_PointsShape.shp",cluster_tolerance="1 Meter") ap.Snap_edit("RoadProximityCBSUsers_PointsShape.shp",[[ap.env.workspace+"//"+RoadNetwork, "EDGE", "0.5 Meter"]]) ap.Snap_edit("RoadProximityCBSUsers_PointsShape.shp",[[ap.env.workspace+"//"+RoadNetwork, "END", "0.5 Meter"]]) ap.AddXY_management("RoadProximityCBSUsers_PointsShape.shp") ap.Dissolve_management(in_features="RoadProximityCBSUsers_PointsShape.shp",out_feature_class="RoadProximityCBSUsers_PointsShape_ConflictFree.shp", dissolve_field=["POINT_X", "POINT_Y"], multi_part="SINGLE_PART") #Set Entrance/Exit Points ap.Buffer_analysis(ExitPoint, "ExitPointA_Buffer.shp","0.2 Meter") ap.Erase_analysis(in_features=RoadVerticesTransferStations, erase_features="ExitPointA_Buffer.shp", out_feature_class="RoadVerticesTransferStations_withoutExit.shp") #Merge all the sets of points together NearPoints=[] ap.Merge_management([ExitPoint,"RoadProximityCBSUsers_PointsShape_ConflictFree.shp", "RoadVerticesTransferStations_withoutExit.shp"], "Final_Points.shp") ap.AddXY_management("Final_Points.shp") NearPoints1=ap.GetCount_management(ExitPoint) NearPoints.append(int(NearPoints1[0])) NearPoints1=ap.GetCount_management("RoadProximityCBSUsers_PointsShape_ConflictFree.shp") NearPoints.append(int(NearPoints1[0])) NearPoints.append(int(RoadVerticesTransferStations_RV_Num)) NearPoints.append(int(RoadVerticesTransferStations_TS_Num)) df = pd.DataFrame(data={"DisMatRows": NearPoints}) df.to_csv(FilePath+"\\ConstructedCSVFiles\\"+Location+"\\"+"CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i) + "_DistanceMatrixComponents.csv", sep=',',index=False) del NearPoints1 ### #Now use these points to split the dissolved road network (defined at start of script) to define the edges ap.SplitLineAtPoint_management(RoadNetwork_Full, "Final_Points.shp", "Final_Lines.shp","0.2 Meters") ap.AddField_management("Final_Lines.shp", "Length", "FLOAT") ap.CalculateGeometryAttributes_management("Final_Lines.shp", [["Length","LENGTH"]],"METERS") ap.GenerateNearTable_analysis("Final_Lines.shp", "Final_Points.shp", "FinalProximityTable.csv", search_radius="0.2 Meter", location="LOCATION", angle="ANGLE", closest="ALL", closest_count=2) #Get walking distances #ap.GenerateNearTable_analysis("CBSUsers_Projected.shp","RoadNetwork_Projected.shp","WalkingDistances_NearTable.csv",location="LOCATION",closest = 'CLOSEST') #Final copying ap.TableToTable_conversion("Final_Lines.shp",FilePath+"\\ConstructedCSVFiles\\"+Location+"\\","CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i) + "_LinesCSV.csv") ap.TableToTable_conversion("Final_Points.shp", FilePath+"\\ConstructedCSVFiles\\"+Location+"\\", "CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i)+ "_PointsCSV.csv") ap.TableToTable_conversion("FinalProximityTable.csv", FilePath+"\\ConstructedCSVFiles\\"+Location+"\\", "CBSUsers_Num_"+str(NumUsers)+"_Set_"+str(i) + "_ProximityTableCSV.csv") bar.update(i) #%% Export Transfer Station Locations #This output is required in the later script for Strategy C ap.AddXY_management(TransferStations) ap.TableToTable_conversion(TransferStations, GeometricOutput_Folder+"\\", "TransferStationsLocations.csv") #################### #%% Export CBS User Locations (from Point) # Exports the nearest points as a shapefile CBSPointsFiles=[] GeoDirectory = FilePath+"\\GeometricOutput\\"+Location for file in os.listdir(GeoDirectory): filename = os.fsdecode(file) if filename.endswith(".shp") and filename.startswith("CBSUsers_Num") : # print(os.path.join(directory, filename)) CBSPointsFiles.append(filename) continue else: continue del file, filename for i in range(0,len(CBSPointsFiles)): Outputtext=CBSPointsFiles[i].replace(".shp","") ap.AddXY_management(FilePath+"\\GeometricOutput\\"+Location+"\\"+CBSPointsFiles[i]) ap.TableToTable_conversion(FilePath+"\\GeometricOutput\\"+Location+"\\"+CBSPointsFiles[i],FilePath+"\\ConstructedCSVFiles\\"+Location,Outputtext+"_CBSPoints.csv")