0

I am trying to code a tool that will evaluate all vertices along a polyline and then find the closest point from a survey shapefile and apply that elevation to that vertex, I am wanting to iterate through a bag of vertices and then at each vertex iterate through the points shapefile and use the distance between the points and elevation in a dictionary and then find the smallest distance at each vertex and grab that elevation here is the code I have so far.

Here is a map package with a typical map https://drive.google.com/file/d/12jrnIn6BAaduLx10_6sI2OJhWB7-IBr2/view?usp=sharing

import os
import arcpy
from math import sqrt
##Declare map and layers, split layer from all layers and then parse .shp    files
##declare map
mxd=arcpy.mapping.MapDocument("CURRENT")
##dataframe
df=arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
 ##pipelayer
listLayers=arcpy.mapping.ListLayers(mxd,"PL_PIPELINE_LN",df)
staging="c:\\staging"
fileName="Report.txt"
path=os.path.join(staging,fileName)
f=open(path,"w+")
##point features to get elevation from
where=""""CODE" = 'BEND' OR "CODE" = 'FLN' OR "CODE" = 'NIP' OR "CODE" = 'RED' OR "CODE" = 'TAP' OR "CODE" = 'TEE' OR "CODE" = 'VLV' OR "CODE" = 'WELD'"""

#lists and dicts for values
pointDict={}
pointList=[]
pointPair=[]
featureList=["FLN","NIP","RED","TAP","TEE","VLV","WELD"]
featurePoints={}
lines=[]

#spatial reference to measure stateplane values
spatialRef=sr_state=arcpy.SpatialReference(6578)

#get centerline name so we can select those features
cenName=""
with arcpy.da.SearchCursor("PL_PIPELINE_LN","CENTERLINE_NAME") as sc:
    for row in sc:
        cenName = row[0]





##get count of features for iteration  , haven't used this yet but may to control iteration
arcpy.SelectLayerByAttribute_management("PL_PIPELINE_LN","NEW_SELECTION","""CENTERLINE_NAME = '{}'""".format(cenName))
result=arcpy.GetCount_management("PL_PIPELINE_LN")
count=int(result.getOutput(0))
print(count)

##clear selection and begin process
arcpy.SelectLayerByAttribute_management("PL_PIPELINE_LN","CLEAR_SELECTION")


##select line features with cen name
with arcpy.da.SearchCursor("PL_PIPELINE_LN",["SHAPE@","OBJECTID"],spatial_reference=spatialRef) as sc:
    for row in sc:
        ##access geometry and get into vertex
        geometry=row[0]
        objId=row[1]
        print (objId)
        pts=geometry.getPart(0)
        for point in pts:
            ##for vertex in vertices
            #print (pointList)
            linePoint=arcpy.PointGeometry(point)
            #featurePoints={}
            ##clear pointDict for each vertex
            pointDict={}
            ##iterate through the survey shape file store values to point dictionary
            with arcpy.da.SearchCursor("1708031_20180611",["SHAPE@","CODE","ELEV","FID"],where,spatial_reference=spatialRef) as sc2:
                for row2 in sc2:
                    pointGeometry= row2[0]
                    code=row2[1]
                    elev=row2[2]
                    oid=row2[3]
                    shpPoint=pointGeometry.getPart(0)
                    ##get distance to vertex
                    dist=sqrt((point.X-shpPoint.X)**2+(point.Y-shpPoint.Y)**2)
                    ##print(dist)
                    ##add values to point dictionary
                    pointDict[oid]=dist,elev,code

            ##find smallest distance in dictioanry
            minDist=(min(pointDict.items(), key=lambda x:x[1]))
            print(minDist[1][0])

        ##minDist=(min(featurePoints.items(),key=lambda x:x[1]))

I get good returns but I don't always get a minimum I suspect my dictionaries aren't complete.

  • I am wondering if I get a point count from the survey shape file and iterate through that it may work better, also I noticed when I traced the distances back to source that I was getting distances from points that shouldn't have been available because of the where statement I included in the SearchCursor, – Gary Lester Jun 20 '18 at 11:04
1

I had 2 problems, one my data set wasn't complete which led to some points along the polyline being very far from the closest survey point

And 2nd problem I needed to take into account scientific notation when dealing with a very small distance. Here is the finished code.

here is working code.

# Name:        module1
# Purpose:
#
# Author:      UseHerName
#
# Created:     17/06/2018
# Copyright:   (c) UseHerName 2018
# Licence:     <your licence>
#-------------------------------------------------------------------------------

import xlwt
from xlwt import easyxf
import os
import arcpy
from math import sqrt

def allPointsLayer(m,d):
    layerListComplete=arcpy.mapping.ListLayers(m,"*",d)
    for layer in layerListComplete:
        if "merged" in layer.name:
            if layer.dataSource.endswith(".shp"):
                allPointsLayerName=layer.name
                return allPointsLayerName


try:
    mxd=arcpy.mapping.MapDocument("CURRENT")
    df=arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
    listLayers=arcpy.mapping.ListLayers(mxd,"PL_PIPELINE_LN",df)[0]
    red_with_white_font = easyxf("pattern: pattern solid, fore_color red; font: color-index white, bold on;")
    allPointsLayer(mxd,df)

except:
    arcpy.AddMessage("You need a pipeline layer or a merged layer")

projectNumber="1708031"
staging="c:\\staging"
fileName="Report.xls"

if os.path.exists(os.path.join(staging,projectNumber)):
    path=os.path.join(staging,projectNumber,fileName)

else:
    os.mkdir(os.path.join(staging,projectNumber))
    path=os.path.join(staging,projectNumber,fileName)

distanceReport=xlwt.Workbook()
worksheet=distanceReport.add_sheet(projectNumber)
worksheet.write(0,0,"Feature OID")
worksheet.write(0,1,"Distance To")
worksheet.write(0,2,"Point Id")
worksheet.write(0,3,"Elevation")
worksheet.write(0,4,"Code")
distanceReport.save(path)
##point features to get elevation from
where=""""CODE" LIKE 'BEND' OR "CODE" LIKE 'CLOSURE' OR "CODE" LIKE 'FLN' OR "CODE" LIKE 'NIP' OR "CODE" LIKE 'RED' OR "CODE" LIKE 'TAP' OR "CODE" LIKE 'TEE' OR "CODE" LIKE 'VLV' OR "CODE" LIKE 'WELD' OR "CODE" LIKE 'VALVE' OR "CODE" LIKE 'ELBOW' OR "CODE" LIKE 'FLANGE'"""
#lists and dicts for values
pointDict={}
pointList=[]
pointPair=[]
featureList=["FLN","NIP","RED","TAP","TEE","VLV","WELD"]
featurePoints={}
lines=[]
allThePoints=allPointsLayer(mxd,df)
#spatial reference to measure stateplane values
spatialRef=sr_state=arcpy.SpatialReference(6578)
#get centerline name so we can select those features
def funCenNAme():
    with arcpy.da.SearchCursor("PL_PIPELINE_LN","CENTERLINE_NAME") as sc:
        for row in sc:
            cenName = row[0]
            return cenName

cenName=funCenNAme()


##get count of features for iteration  , haven't used this yet but may to control iteration
arcpy.SelectLayerByAttribute_management("PL_PIPELINE_LN","NEW_SELECTION","""CENTERLINE_NAME = '{}'""".format(cenName))
result=arcpy.GetCount_management("PL_PIPELINE_LN")
count=int(result.getOutput(0))

##clear selection and begin process
arcpy.SelectLayerByAttribute_management("PL_PIPELINE_LN","CLEAR_SELECTION")

i=1
##select line features with cen name
with arcpy.da.UpdateCursor("PL_PIPELINE_LN",["SHAPE@","OBJECTID"],spatial_reference=spatialRef) as sc:
    for row in sc:
        ##access geometry and get into vertex
        geometry=row[0]
        objId=row[1]
        pts=geometry.getPart(0)
        pointArray=arcpy.Array()
        for point in pts:
            i+=1
            sheetRow=worksheet.row(i)
            ##for vertex in vertices
            #print (pointList)
            linePoint=arcpy.PointGeometry(point)
            #featurePoints={}
            ##clear pointDict for each vertex
            pointDict={}
            newPoint=arcpy.Point()
            ##iterate through the survey shape file store values to point dictionary
            with arcpy.da.SearchCursor("1708031_merged",["SHAPE@","CODE","ELEV","FID"],where,spatial_reference=spatialRef) as sc2:
                for row2 in sc2:
                    pointGeometry= row2[0]
                    code=row2[1]
                    elev=row2[2]
                    oid=row2[3]
                    shpPoint=pointGeometry.getPart(0)
                    ##get distance to vertex
                    dist=sqrt((point.X-shpPoint.X)**2+(point.Y-shpPoint.Y)**2)
                    ##print(dist)
                    ##add values to point dictionary
                    pointDict[oid]=dist,elev,code,oid

            print point.X
            print point.Y
            print minDist[1][1]
            minDist=(min(pointDict.items(), key=lambda x:x[1]))
            newPoint.X=point.X
            newPoint.Y=point.Y
            newPoint.Z = minDist[1][1]
            pointArray.add(newPoint)
            if round(minDist[1][0],4)>.300:
                worksheet.write(i,0,objId)
                print ("OID For Feature IS {0} And the vertex nummber is {1}".format(objId,i))
                worksheet.write(i,1,round(minDist[1][0],4),red_with_white_font)
                print ("""FID of elevation point {}""".format(str(minDist[1][3]).encode('utf8')))
                worksheet.write(i,2,minDist[1][3],red_with_white_font)
                print ("""Elevation: {}""".format(str(minDist[1][1]).encode('utf8')))
                worksheet.write(i,3,float(minDist[1][1]),red_with_white_font)
                print ("""Distance to point {}""".format(str(round(minDist[1][0],4)).encode('utf8')))
                print (minDist[1][2])
                worksheet.write(i,4,minDist[1][2],red_with_white_font)
            else:
                worksheet.write(i,0,objId)
                worksheet.write(i,1,round(minDist[1][0],4))
                worksheet.write(i,2,minDist[1][3])
                worksheet.write(i,3,float(minDist[1][1]))
                worksheet.write(i,4,minDist[1][2])

        newGeometry=arcpy.Polyline(pointArray,None,True)
        print newGeometry.centroid.X
        row=newGeometry,objId
        sc.updateRow(row)

distanceReport.save(path)
|improve this answer|||||

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.