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!")

ArcIdeas: Scale-Dependant Settings.

Someone mentioned an idea on ArcIdeas for making various display settings on a feature classes scale-dependent.  Right now some of that can be accomplished by loading a feature classes multiple times, adjusting the settings, and setting the visible range.  Working more and more in ArcGIS Server, I can see the value of increased scale-dependent settings.

I’m not sure how rapidly ESRI takes “Ideas” into consideration but if you feel like it would benefit you, why not promote this idea:  Scale Range, SQL Query and Symbology Rendering in ArcMap.

Have a good weekend, all!

Using arcpy to List Domains Assigned to Featureclass Fields

I was making an edit (adding leading “0”s) to a coded-value domain in an SDE database and realized that my edits were changing the order of the rows of my domain.  Rows were moved to the bottom of the list when they were edited.

So I went through the process of converting my domain back to a table, made my edits in Access and exported the rows to a .dbf in the order I wanted them.

I removed the domain from the field I knew it was assigned to but when I tried to delete the domain, I received an error (The domain is used as a default domain).

The Google Machine led me to an ArcForums post by  with some python code for listing all the fields with a domain.

I tweaked the original a bit, first because it was only examining feature classes in a feature dataset, not stand-alone feature datasets.  And secondly, I updated it to use arcpy.  I posted both the 9.3 version and the 10.0 version, but here is a quick look.  The key addition is the ‘listfc(“”)’ line that is the first line of the def listds() module.

import arcpy

min_workspace = r"C:\Users\mranter\AppData\Roaming\ESRI\Desktop10.0\ArcCatalog\min.minstaff.sde"

infgdb=(min_workspace)
arcpy.env.workspace = infgdb

def listfc(inDataset):
   featureclasses = arcpy.ListFeatureClasses("","",inDataset)
   for f in featureclasses:
      print "feature class: ",f

lfields=arcpy.ListFields(f)

for lf in lfields:
   if lf.domain < "":
      print "      domain",f, lf.name, lf.domain

def listds():
   listfc("")

   datasets=arcpy.ListDatasets ("","")
   for d in datasets:
      print "  dataset: ",d

listfc(d)
listds()

Revoking ‘Public’ Permissions in PostgreSQL ArcSDE

While banging my head on how to grant access to a referenced mosaic dataset , I did something out of frustration that I normally would not do–I granted ‘public’ access to some data.

Then, after figuring out the problem, I went to revoke public access using ArcCatalog and received this error message:

Error 999999: Error executing function.

The Object being referenced does not exist [ERROR:  role “public” does not exist::SQL state: 42704]
Failed to execute (ChangePrivileges).

I actually expected that would occur–ArcSDE isn’t aware of PostgreSQL’s public account.

The solution was to go into pgAdmin II and revoke the privileges there.

Extract Values to Points (Spatial Analyst) Bug

One of the Spatial Analyst tools we often use in ArcGIS is the “Extract Values to Points” tool.  This allows us to take a point file (well locations in our case) and attach a value (elevations) from a raster image (a DEM) to each point.

Today I was running it for the first time against an Image Service we recently published and I received a warning message,”WARNING 000957: Skipping feature(s) because of NULL or EMPTY geometry”.  But the script seemed to run and the final results said “Succeeded” so I thought it was probably fine.  But as I double-checked, I realized the results were wonky.

Turns out that I had two records with Null geometry in my point file of 397 records.  These two records threw the above error but actually had a value in the [Rastervalu] field.  Turns out all 397 records had values.  These two records were consecutive–let’s say the 100th and 101st records in my shapefile.  What happened is record 100 got the value for record 102, record 101 get the value for 103, record 102 (which has valid geometry) had the value for 104.  This pattern, each record having the value for the record 2 place after it, continued until record 396 which had the value for record 397.  Record 397 also had the value for 397.  So the final three records all had the value for the final record.

What I would have expected would be for the two records with Null Geometry to have null values in the [Rastervalu] field and the rest of the records to have the correct values.  Despite the warning, it is very misleading for all the records to end up with a value.

I have a simplified example below.  I made a point shapefile with four records.  The first, third, and fourth  records have valid geometries; the second has Null geometry.  The second record ends up with the value for the third record.  The third record, has the value for the fourth.  The fourth record being the last record, ends up with the last valid value, which was its own.

The results that I would have hoped for would be for the third record to have a Null value.

The way I envision what is occurring behind the scenes is this:  the process makes a list (more of a stack in programming terms) of result values as it processed the points but just assumes that every record will return a value so it does not track which value goes with which shape.

When it reached the two null geometries, it threw an error but continued on.  It did not add a value for these records to the stack of values–when it comes across records with valid geometry but do not intersect the raster it adds a psuedo-null value of (-9999) to the stack.  After it processed all the records it had 395 values in the stack.  It then went, one-by-one through the stack and populated the records in the output shapefile, the first record got the first value in the stack, the second record got the second value, the 100th record got the 100th value (which came from the location of the 102nd record) and so on.  At the end, the final two records received the last valid value.

This final behavior–using the last valid value–corresponds a bit to a behavior we’ve seen with ArcObjects in general.  When iterating through a table, if a field is Null for a specific row, the value from the last non-Null value for that field is often returned.

I’m in the process of submitting a bug to ESRI.  I’m not sure if this existed prior to ArcGIS 10.0 (I’m guessing it did) or if it occurs in other processes (I’m guessing it does).  I did find out that the “Extract Multi Values to Points” works as expected.  I’m guessing it is because unlike the “Extract Values to Points” which creates a new shapefile, this tool appends fields to the existing shapefile and presumably processes records one-by-one without putting the results in a virtual stack.  The “Extract Multi Values to Points” tool also does not throw any warnings.


ArcGIS 9.3.1 to ArcSDE 10 Connection Error

We finally installed an instance of ArcSDE 10 today.  My first attempt at connecting in ArcCatalog 9.3.1 failed with the following error:

Failed to connect to the specified server.

This release of the GeoDatabase is either invalid or out of date. [Please run the

ArcSDE setup utility using the -o install option.]

DBMS table not found[sde.sde.GDB_Release]

 

Turns out the solution was simple, this article points out that Service pack 2 is required.  After updating my ArcGIS, the problem was solved:

I was the proud, new recipient of a non-misleading error message.  After some minor head-slapping, I realized that 9.3.1 and older versions of ArcGIS are not able to connect to ArcSDE 10.  Doh!  My bad–I should have read ESRI’s clear information on this topic.  At least I never actually tried to install ArcSDe with the -o option as the first error message instructed me to do.

Launching a Python script with parameters–3 methods.

Since I use python for different tasks, I launch python scripts a variety of ways. Depending on what I am doing, a single script may need to accept parameters from either:

  1. Passed in from an ArcGIS Toolbox Tool.
  2. Re-occurring default value.  Often used in scheduled processes, a nightly backup, for example.
  3. A temporary set of values used in an interactive, debugging session.

What I often do is make the parameter interpretation flexible to meet my needs.  The sample below shows how I do this.  The logic first checks to see if the correct number of parameters were used to launch the script (i.e. if it is called by an Arc Toolbox tool), where the files for the default files exist or if the debug values are valid, including checking the current date against a hard-coded date variable.  Juggling the conditional structure would allow you to prioritize the options differently.

I am also using Tkinter to display an interactive dialog if none of the three conditions are successfully met.

import os.path
import datetime
import shutil
import sys

theEmailText = "nStart: " + str(datetime.datetime.now())

#This example shows three different ways for a script to receive input paramters.
#
# First, if it was run with the necessary number of parameters, as if launched by
# an ArcToolbox tool, it uses those parameters.

if len(sys.argv) == 2:
    wellsShapeFile = sys.argv[0]
    unwellsShapeFile = sys.argv[1]
    theEmailText = theEmailText + "nUsing Parameters Passed In:nn Located Wells File: "+wellsShapeFile+"n Unlocated Wells File: "+unwellsShapeFile
else:
    dateString = datetime.date.today().strftime("%Y%m%d")
    wellsShapeFile = "C:/cwi5_bk/wells/temp/wells_.shp"
    unwellsShapeFile = "C:/cwi5_bk/wells/temp/unloc_wells.shp"

    #Second attempt is if there are default values that should be used.
    #I use this for a process that is run via scheduled Windows Task
    if (os.path.exists(wellsShapeFile)) and (os.path.exists(unwellsShapeFile)):
        theEmailText = theEmailText + "nUsing Automated Date-Based file names:nn Located Wells File: "+wellsShapeFile+"n Unlocated Wells File: "+unwellsShapeFile
    else:
        #The third method is used for debugging/running in the IDE.
        #I put a check condition so this is valid for one day only
        #And then hard-code the temporary paths.
        #
        #Note that you may want to modify the structure of the IF statement
        #used for methods 2 & 3 so that it checks for the manual-override (3rd)
        #method first
        if dateString == '20101213':
            wellsShapeFile = "C:/cwi5_bk/wells/temp/wells.shp"
            unwellsShapeFile = "C:/cwi5_bk/wells/temp/unloc.shp"
            theEmailText = theEmailText + "nUsing Manual Override file names:nn Located Wells File: "+wellsShapeFile+"n Unlocated Wells File: "+unwellsShapeFile
        else:
            theEmailText = theEmailText + "n Manual Override for file names does not meet data filter"

    if (not os.path.exists(wellsShapeFile)) or (not os.path.exists(unwellsShapeFile)):
        from Tkinter import *
        msgbox = Tk()
        msgbox.title('Error')
        Message(msgbox,text="Must either use ArcTool to launch or edit file parameters", bg='royalblue',fg='ivory', relief=GROOVE).pack(padx=10, pady=10)
        msgbox.mainloop()
        print theEmailText
        quit()

print theEmailText

Another TopoToRaster Error

Subtitled:  Why error messages are good.

Came up with another error while running TopoToRaster but this time ArcGIS gave an error message that led to a solution.  Turned out all my contour lines had an elevation of 16 which TopoToRaster did not like.  I had intended to increase the elevation and inadvertently set them all to sixteen.  I had saved the previous values before editing so it turned out to be a simple fix and I didn’t have to spend a day trying figure out what was wrong.


Trying to perform “relate” in Model Builder

I was using Model Builder (ugh!) to select records in one table (CWI.C5ST) that relate to a subset of records ([BHGEOPHYS] = ‘Y’) in another table (CWI.C5IX).  There is not an existing tool for doing this in ArcGIS.  I did find a post by Layne Seely in ArcForums titled “trying to perform “relate” in Model Builder.” that led me to the Make Table View under Data Management Tools-Make Table View and even had  the basic syntax I was looking for (if it hadn’t I probably would have guessed that it would not allow subqueries).

Since I always like to see exactly the syntax someone uses, I decided to show the syntax I used in case that helps anyone else.

Sample of a Table View using a subquery

It does require that the tables reside in the same workspace, which is a reasonable limitation that I can live with.  Another drawback is that if you export the model to a python script, the command can get pretty long, especially if you have many fields.

ArcGIS Ideas

Ok, maybe I have been hiding under a rock, but just saw this is a post from James Fee’s Spatially Adjusted blog.  ESRI has a public suggestion box called  ArcGIS Ideas where users can submit suggestions for ESRI.  One thing I like is that it is open for other users to see and comment.  Not sure how responsive ESRI will be but at least it gives us a place to vent.