Quick & Dirty arcpy: Batch Splitting Polylines to a Specific Length.

For some odd reason, I wanted to split all the arcs in a polyline feature class to a specific length–if a specific feature was longer than the target length, it would become two or more separate polyline records.

Here is the bare-bones script that copies an existing feature class into a new feature class then processes each record, splitting it into multiple records if the polyline is longer than the user-specified tolerance.  Some cautionary notes:

  • This is Quick & Dirty code–minimal error catching or documentation.
  • I basically tested this against one feature class (the one I wanted to split) once I got it to work, I quit.
  • There is some rounding error–features may be a tad bit off (a few ten-thousandths of a unit).
  • I did not test against multi-part features.
  • The tolerance is the native units of the data–if your data is in meters but you want to split the polylines every mile, enter 1,609.344.

I have included both a toolbox file (.tbx) and python script (.py).  After loading the toolbox, you’ll have to change the Source of the script by right-clicking on it, selecting the Source tab, and then navigating to the .py file.

Here is the code for the Googlebots, but you are better off just downloading it.

import arcpy
import sys, math

def printit(inMessage):
    print inMessage
    arcpy.AddMessage(inMessage)

if len(sys.argv) > 1:
    inFC = sys.argv[1]
    outFC = sys.argv[2]
    alongDistin = sys.argv[3]
    alongDist = float(alongDistin)
else:
    inFC = "C:/temp/asdfasdf.mdb/jkl"
    OutDir = "C:/temp/asdfasdf.mdb"
    outFCName = "jkl2d"
    outFC = OutDir+"/"+outFCName
    alongDist = 1000

if (arcpy.Exists(inFC)):
    print(inFC+" does exist")
else:
    print("Cancelling, "+inFC+" does not exist")
    sys.exit(0)

def distPoint(p1, p2):
    calc1 = p1.X - p2.X
    calc2 = p1.Y - p2.Y

    return math.sqrt((calc1**2)+(calc2**2))

def midpoint(prevpoint,nextpoint,targetDist,totalDist):
    newX = prevpoint.X + ((nextpoint.X - prevpoint.X) * (targetDist/totalDist))
    newY = prevpoint.Y + ((nextpoint.Y - prevpoint.Y) * (targetDist/totalDist))
    return arcpy.Point(newX, newY)

def splitShape(feat,splitDist):
    # Count the number of points in the current multipart feature
    #
    partcount = feat.partCount
    partnum = 0
    # Enter while loop for each part in the feature (if a singlepart feature
    # this will occur only once)
    #
    lineArray = arcpy.Array()

    while partnum < partcount:
        # Print the part number
        #
        #print "Part " + str(partnum) + ":"
        part = feat.getPart(partnum)
        #print part.count

        totalDist = 0

        pnt = part.next()
        pntcount = 0

        prevpoint = None
        shapelist = []

        # Enter while loop for each vertex
        #
        while pnt:

            if not (prevpoint is None):
                thisDist = distPoint(prevpoint,pnt)
                maxAdditionalDist = splitDist - totalDist

                print thisDist, totalDist, maxAdditionalDist

                if (totalDist+thisDist)> splitDist:
                    while(totalDist+thisDist) > splitDist:
                        maxAdditionalDist = splitDist - totalDist
                        #print thisDist, totalDist, maxAdditionalDist
                        newpoint = midpoint(prevpoint,pnt,maxAdditionalDist,thisDist)
                        lineArray.add(newpoint)
                        shapelist.append(lineArray)

                        lineArray = arcpy.Array()
                        lineArray.add(newpoint)
                        prevpoint = newpoint
                        thisDist = distPoint(prevpoint,pnt)
                        totalDist = 0

                    lineArray.add(pnt)
                    totalDist+=thisDist
                else:
                    totalDist+=thisDist
                    lineArray.add(pnt)
                    #shapelist.append(lineArray)
            else:
                lineArray.add(pnt)
                totalDist = 0

            prevpoint = pnt                
            pntcount += 1

            pnt = part.next()

            # If pnt is null, either the part is finished or there is an
            #   interior ring
            #
            if not pnt:
                pnt = part.next()
                if pnt:
                    print "Interior Ring:"
        partnum += 1

    if (lineArray.count > 1):
        shapelist.append(lineArray)

    return shapelist

if arcpy.Exists(outFC):
    arcpy.Delete_management(outFC)

arcpy.Copy_management(inFC,outFC)

#origDesc = arcpy.Describe(inFC)
#sR = origDesc.spatialReference

#revDesc = arcpy.Describe(outFC)
#revDesc.ShapeFieldName

deleterows = arcpy.UpdateCursor(outFC)
for iDRow in deleterows:       
     deleterows.deleteRow(iDRow)

del iDRow
del deleterows

inputRows = arcpy.SearchCursor(inFC)
outputRows = arcpy.InsertCursor(outFC)
fields = arcpy.ListFields(inFC)

numRecords = int(arcpy.GetCount_management(inFC).getOutput(0))
OnePercentThreshold = numRecords // 100

printit(numRecords)

iCounter = 0
iCounter2 = 0

for iInRow in inputRows:
    inGeom = iInRow.shape
    iCounter+=1
    iCounter2+=1    
    if (iCounter2 > (OnePercentThreshold+0)):
        printit("Processing Record "+str(iCounter) + " of "+ str(numRecords))
        iCounter2=0

    if (inGeom.length > alongDist):
        shapeList = splitShape(iInRow.shape,alongDist)

        for itmp in shapeList:
            newRow = outputRows.newRow()
            for ifield in fields:
                if (ifield.editable):
                    newRow.setValue(ifield.name,iInRow.getValue(ifield.name))
            newRow.shape = itmp
            outputRows.insertRow(newRow)
    else:
        outputRows.insertRow(iInRow)

del inputRows
del outputRows

printit("Done!")

Quick & Dirty arcpy: Field Listings

I have to often get a table structure for a feature class or table into either a spreadsheet or word processing document.  There might be an easy way to do this in ArcGIS 10 but I haven’t found it.  So, as is my nature, I decided to roll my own.

This is a bare-bones script that iterates through the fields, printing the field name, type, width, and precision.  There are three optional features to it:

  • You can choose to have it list the domain, if there is one, on each field.
  • You can have it write to a text file (otherwise you can just copy & paste the results from the results window).
  • You can have it count the number of populated records.  This can take a long time if working with a large dataset.  Also note that my logic for determining what constitutes being populated may not be what you need but the structure is there.  I also do not account for all field types, if the field is of a type I have not account for, the code will return -999.

To use the script from ArcToolbox, you need to pass it four parameters, their Names, type, whether they are input or output, and whether they are required or optional are:

  • featureclass, Table, Input, Required
  • includedomainstring, Boolean, Input, Required (controls whether or not domains are exported)
  • doCountsRespone, Boolean, Input, Required (controls whether or not you want to get the number of populated records.  (Your definition of populated may vary from my code)
  • outputFile, File, Output, Optional (optional output file to write)

Here is the code, but you are better off just downloading it since I haven’t figured out a good way to have WordPress play nicely with python’s indenting.

# Name: ListFields-arcpy.py
#
# Purpose: Lists the fields, their type, width, and precision
# Can either have it export it to a CSV file or copy
# and paste from the results window.
#
# To use, create a tool from the script and add 3 parameters:
#  1) Table, Input, Required
#  2) Boolean, Input, Required (controls whether or not domains are exported)
#  3) Boolean, Input, Rekquired (controls whether or not you want to get the number of
#  Populated records.&nbsp; (Your defintion of populated may vary from my code)
#  4) File, Output, Optional (optional output file to write)
#
#

import arcpy,sys,os

def printit(inMessage):
 print inMessage
 arcpy.AddMessage(inMessage)

if len(sys.argv) > 4:
 featureclass = sys.argv[1]
 includedomainstring = sys.argv[2]
 doCountsRespone = sys.argv[3]
 outputFile = sys.argv[4]
else:
 featureclass = "C:/temp/before.shp"
 includedomainstring = "false"
 doCountsRespone = "true"
 outputFile = "C:/temp/before.csv"

if (outputFile == ""):
 doOutputFile = False
else:
 doOutputFile = True

if (str(doCountsRespone).lower() == "true"):
 doCounts = True
else:
 doCounts = False

if (str(includedomainstring).lower() == "true"):
 includedomain = True
else:
 includedomain = False

lfields=arcpy.ListFields(featureclass)

d = arcpy.Describe(featureclass)
printit("Dataset: "+d.baseName)
printit("Type: "+d.dataType)
printit("Path: "+d.catalogPath)
printit(" ")

tableheaders = 'name,type,width,precision'

if (doCounts == True):
 tableheaders+=",count"

if (includedomain == True):
 tableheaders+=",domain"

if (doOutputFile):
 tmpfile = open(outputFile,"w")
 tmpfile.write(tableheaders)
 tmpfile.write("n")

printit (tableheaders)
for lf in lfields:

 pThisline = lf.name+","+lf.type +","+str(lf.length)+","+str(lf.precision)

 if (doCounts == True):

 rowCount = 0

 #Note that I do not account for all field types
 #Also note that my definition of being populated may vary from yours.
 #I am using -999 as a flag to indicate a field type was not successfully
 #identified.
 if (lf.type == "Double") or (lf.type == "Single")&nbsp; or (lf.type == "Integer") or (lf.type == "SmallInteger"):
  queryString = '"'+lf.name + '" > 0'
  rows = arcpy.SearchCursor(featureclass, queryString, "", "", "")
 elif (lf.type == "String"):
  queryString = '"'+lf.name + '" <> ' + "''"
  rows = arcpy.SearchCursor(featureclass, queryString, "", "", "")
 else:
  rowCount = -999
  #rows = arcpy.SearchCursor(featureclass, "", "", "", "")

 if (rowCount == 0):
  for row in rows:
   rowCount+=1

 pThisline=pThisline+","+str(rowCount)

 if (includedomain == True):
  pThisline=pThisline+","+lf.domain

 printit (pThisline)

 if (doOutputFile):
  tmpfile.write(pThisline)
  tmpfile.write("n")

if (doOutputFile):
 tmpfile.close

Multiple outputs for Python scripts

Related to my post on how I enable a script to accept parameters from different sources, I also often set up pythons scripts to output information a variety of ways.  This is largely due to the fact that some are called by ArcToolbox scripts.  Running in ESRI’s domain, these scripts need to send the output through the arcgisscripting object but if you are running the python outside the ArcGIS framework, you can just print.

If you assume one output method but then run your code in the opposite framework, you don’t get to see all the pretty little messages.  What I do is create a simple little routine that broadcasts the message both ways.  This is probably an obvious solution but took a few cases before I went ahead and started implementing it.

gp = arcgisscripting.create()

#This will print both to the geoprocessing window or Python output window
def gpprint(inmessage):
 gp.addmessage(inmessage)
 print inmessage

<Code to do stuff>

#Ok, I want to send a message:
gpprint('Hello, sailor!')

Zipping a shapefile via ArcToolbox

UPDATE:

After receiving a request to modify the code to ignore .lock files, I have an updated to this post.

 

I’ve received a request on how to use the Zip Shapefile code I posted last week from ArcGIS. Sorry, I did not set the code up to call directly from ArcGIS but only as an illustration of how it can be done.

I have, however, with some minor tweaking, made a version that can added to ArcToolbox. The steps to install are below, please note that at one point I had you download a *.zip file that had been renamed to *.jpg–this should now be corrected and the link should lead you directly to zipshapefile.zip.  Because of this steps two and three are obsolete.

1) Download the code from here.
2) Rename the file from zipshapefile-zip.jpg back to zipshapefile.zip.
3) Unzip the file.
4) Move ZipShapefile.py to C:Program FilesArcGISArcToolBoxScriptsZipShapefile.py.
5) Optionally, move Zip Shapefile.tbx, perhaps C:Program FilesArcGISArcToolBoxToolboxes.
6) Add the toolbox to ArcToolbox. ESRI has instructions here on how to do this.

You should now have a new toolbox named “Zip Shapefile” with a script named “Zip a Shapefile” in it. Clicking on on the tool will bring up this dialog.

**********************************
In response to Chris:

I believe you need to copy the ZipShapefile.py file from the .zip that you downloaded to C:Program FilesArcGISArcToolBoxScripts, the error message is consistent with the tool not being about to find the python script there.

If you prefer to place the ZipShapefile.py in a different location, you will need to change the source on the tool. To do this, right click on the tool in ArcCatalog and change the path of the Script File as set in the Source tab (see below):

NED Processing in Python & Geoprocessing

The USGS updated a significant portion of NED data for my state in June.  I recently downloaded the updates, processed them–projecting and converting the elevations from meters to feet–using python and geoprocessing.  My python skills are still pretty crude but I was able to get the job done.

One of the benefits of working for the public sector is that I can more freely publish code without worrying about “trade secrets” or what have you so I thought I would put my code out for anyone to see, maybe someone will find it useful.  It is split into two separate files, mostly because I used Model Builder to generate the main processing chunk of code (Project_Reclassify.py) and another script to control looping, etc. The code isn’t pretty and you’ll see the results of my development process with two subroutines in LaunchScript.py that could really be one.

Continue reading “NED Processing in Python & Geoprocessing”