Monday, April 15, 2019

ArcGIS Split Line at Vertices without an Advanced License

ArcGIS's Split Line At Vertices is a handy tool, but it requires an advanced license. Luckily with a bit of Python magic you can create a feature class of lines split at vertices just like the tool would create.

The code:


##split line at point script
##Emil Brundage
##emilharold@gmail.com

### USER INPUTS

#input feature class
inFc = r"C:\Path\To\Data.gdb\InPolylines"
#output feature class
#will be created by script
outFc = r"C:\Path\To\Data.gdb\OutPolylines"

### END USER INPUTS

#----------

print "importing"
import arcpy
import os

print "creating new feature class"
outPath, outName = os.path.split (outFc)
outFc = arcpy.CreateUniqueName (outName, outPath)
outPath, outName = os.path.split (outFc)
outFc = arcpy.CreateFeatureclass_management (outPath, outName, "POLYLINE", inFc,
                                             spatial_reference = inFc) [0]

print "determining fields"
flds = [f.name for f in arcpy.ListFields (inFc)]
flds = [f.name for f in arcpy.ListFields (outFc) if f.name in flds] + ["SHAPE@"]

print "iterating"
with arcpy.da.InsertCursor (outFc, flds) as iCurs:
    with arcpy.da.SearchCursor (inFc, flds) as sCurs:
        for row in sCurs:
            row = list (row)
            geom = row [-1]
            if not geom:
                iCurs.insertRow (row)
                continue
            for i in range (geom.partCount):
                part = geom.getPart (i)
                fPoint = None
                count = part.count
                for c in range (count):
                    if fPoint == None:
                        fPoint = part.getObject (c)
                        continue
                    lPoint = part.getObject (c)
                    ln = arcpy.Polyline (arcpy.Array ([fPoint, lPoint]))
                    row [-1] = ln
                    iCurs.insertRow (row)
                    fPoint = lPoint

print
print "created:"
print outFc
print
print "done"
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Tuesday, August 15, 2017

ArcGIS Point Distance Without an Advanced License

The point distance tool generates a table that contains data on distances between points of two feature classes. Below is a function that will create the same table and doesn't require an advanced ArcGIS license.
The code:


import os
import arcpy

def PointDistance (inFeat, nearFeat, outTable):
    #create table
    tabPath, tabName = os.path.split (outTable)
    arcpy.CreateTable_management (tabPath, tabName)
    #add fields
    fldDi = {"INPUT_FID" : "LONG",
             "NEAR_FID" : "LONG",
             "DISTANCE" : "DOUBLE"}
    for fld in ["INPUT_FID", "NEAR_FID", "DISTANCE"]:
        arcpy.AddField_management (outTable, fld, fldDi [fld])
    with arcpy.da.InsertCursor (outTable, ["INPUT_FID",
                                           "NEAR_FID",
                                           "DISTANCE"
                                           ]) as iCurs:
        with arcpy.da.SearchCursor (inFeat,
                                    ["OID@",
                                     "SHAPE@"
                                     ]) as inCurs:
            with arcpy.da.SearchCursor (nearFeat,
                                        ["OID@",
                                         "SHAPE@"]
                                        ) as nearCurs:
                for inOid, inGeom in inCurs:
                    for nearOid, nearGeom in nearCurs:
                        row = (inOid, nearOid,
                               inGeom.angleAndDistanceTo (nearGeom) [1])
                        iCurs.insertRow (row)
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Tuesday, April 25, 2017

Using Field Info to Limit Fields When Field Mapping Isn't Available

I work with feature classes that can have several dozens of fields, sometimes numbering in the triple digits. Once I start a geoprocessing workflow, I'm often only interested in one or two fields from my input feature classes for my analysis. Usually field mappings are the way to go to limit fields. However, certain tools don't have field mappings as an input, so on its own the user is stuck keeping all fields of all inputs involved. Instead of deleting unwanted fields after a geoprocessing it's possible to limit fields by first creating layers with hidden fields using field info. In the below example I've written a function that creates feature layers with unwanted fields set to hidden. I then use these layers in an intersect, which produces a feature class whose attribute table contains only the fields that were set to visible.
The code:

import arcpy

def ApplyFldInfo (fc, flds):
    """
    create layer with desired fields hidden
    fc: the full path of the input feature class
    flds: a python list of fields to remain visible
    emilharold@gmail.com
    """
    #create layer with unique name
    layer = "layer"
    i = 0
    while arcpy.Exists (layer):
            layer = "layer_{}".format (i)
            i += 1
    arcpy.MakeFeatureLayer_management (fc, layer)
    #get field info
    fldInfo = arcpy.Describe (layer).fieldInfo
    #set fields no in input list to hidden
    for i in range (fldInfo.count):
            fldName = fldInfo.getFieldName (i)
            if not fldName in flds:
                    fldInfo.setVisible (i, "HIDDEN")
    #create new layer with field info applied
    outLayer = "outLayer"
    i = 0
    while arcpy.Exists (outLayer):
            outLayer = "outLayer_{}".format (i)
            i += 1
    arcpy.MakeFeatureLayer_management (layer, outLayer,
                                       field_info = fldInfo)
    #delete initial layer
    arcpy.Delete_management (layer)
    return outLayer

#first feature class
fc1 = r"C:\Path\To\Fc1"
#keep field in first feature class
fld1 = "field1"
#second feature class
fc2 = r"C:\Path\To\Fc2"
#keep fields in second feature class
fld2 = "field2"
fld3 = "field3"

#create layers with fieldinfo
layer1 = ApplyFldInfo (fc1, [fld1])
layer2 = ApplyFldInfo (fc2, [fld2, fld3])

#intersect
outFc = r"C:\Path\To\Output"
arcpy.Intersect_analysis ([layer1, layer2], outFc)

#clean up
for layer in [layer1, layer2]:
    arcpy.Delete_management (layer)
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Friday, April 21, 2017

A More Robust Merge

Merging various feature classes into a single feature class is a common GIS task, and any GIS newbie knows that ArcGIS has a Merge tool to get the job done. The tool works well, with one small hiccup. If two features classes being merged have fields of the same name but incompatible types (integer and string, say) the tool will fail. This was probably a conscious decision by the developers of ArcGIS, but I would rather have a merge that makes the best choice for a field type when a conflict arises. In the below function, a list of feature classes is merged and fields are evaluated to determine the proper field type for the output feature class. For example, when two feature classes have the same field name but one is text while the other is numeric, the script creates a text field and converts values from the numeric field to text. Date fields that conflict with other field types are converted to strings. Conflicts between double fields and integer fields default to double. The function's mergeFcs input is a python list of input feature classes. The outFc input is the full path of the output feature class.

The code:

import arcpy
import os

def Merge (mergeFcs, outFc):

    """
    Merge feature classes with field type discrepancies
    emilharold@gmail.com
    """
    
    print "iterating fields to determine field types"
    stringFlds = []
    lIntFlds = []
    floatFlds = []
    doubleFlds = []
    dateFlds = []
    numFcs = len (mergeFcs)
    digits = len (str (numFcs))
    lenDi = {}
    placeDi = {}
    fObDi = {}
    place = 0
    
    for i, fc in enumerate (mergeFcs, 1):
        if i != 1 and not (i - 1) % 10:
            print
        print str (i).zfill (digits) + "/" + str (numFcs),
        for fld in arcpy.ListFields (fc):
            fType = fld.type
            fName = fld.name
            if fType in ("Geometry",
                         "Guid",
                         "GlobalID",
                         "OID") or fld.required:
                continue
            fObDi [fName] = fld
            if not fName in placeDi.values ():
                place += 1
                placeDi [place] = fName
            if fType == "String":
                stringFlds += [fName]
                fLength = fld.length
                if not fName in lenDi:
                    lenDi [fName] = fLength
                elif fLength > lenDi [fName]:
                    lenDi [fName] = fLength
            elif fType == "Integer":
                if fName in dateFlds:
                    stringFlds += [fName]
                    fLength = 50
                    if not fName in lenDi:
                        lenDi [fName] = fLength
                    elif fLength > lenDi [fName]:
                        lenDi [fName] = fLength
                else:  
                    lIntFlds += [fName]
            elif fType == "Single":
                if fName in dateFlds:
                    stringFlds += [fName]
                    fLength = 50
                    if not fName in lenDi:
                        lenDi [fName] = fLength
                    elif fLength > lenDi [fName]:
                        lenDi [fName] = fLength
                else:
                    floatFlds += [fName]
            elif fType == "Double":
                if fName in dateFlds:
                    stringFlds += [fName]
                    fLength = 50
                    if not fName in lenDi:
                        lenDi [fName] = fLength
                    elif fLength > lenDi [fName]:
                        lenDi [fName] = fLength
                else:
                    doubleFlds += [fName]
            elif fType in ("Date", "Blob", "Raster"):
                if fName in lIntFlds + floatFlds + doubleFlds:
                    stringFlds += [fName]
                    fLength = 50
                    if not fName in lenDi:
                        lenDi [fName] = fLength
                    elif fLength > lenDi [fName]:
                        lenDi [fName] = fLength
                else:
                    dateFlds += [fName]
                    

    print
    stringFlds = set (stringFlds)
    lIntFlds = set (lIntFlds)
    floatFlds = set (floatFlds)
    doubleFlds = set (doubleFlds)

    print "creating feature class"
    outLoc, outName = os.path.split (outFc)
    checkFc = mergeFcs [0]
    sr = arcpy.Describe (checkFc).spatialReference
    shapeType = arcpy.Describe (checkFc).shapeType.upper ()
    arcpy.CreateFeatureclass_management (outLoc, outName,
                                         shapeType,
                                         spatial_reference = sr)
    lenFlds = str (len (placeDi.keys ()))
    digits = len (lenFlds)
    print "adding fields"
    for i, place in enumerate ((placeDi.keys ()), 1):
        if i != 1 and not (i - 1) % 5:
            print
        print str (i).zfill (digits) + " of " + lenFlds, "",
        fName = placeDi [place]
        fld = fObDi [fName]
        if fName in stringFlds:
            arcpy.AddField_management (outFc, fName, "TEXT",
                                       field_length =
                                       lenDi [fName])
        elif fName in floatFlds:
            arcpy.AddField_management (outFc, fName, "FLOAT")
        elif fName in doubleFlds:
            arcpy.AddField_management (outFc, fName, "DOUBLE")
        elif fName in lIntFlds:
            arcpy.AddField_management (outFc, fName, "LONG")
        else:
            arcpy.AddField_management (outFc, fName, fld.type)

    print
    print "iterating feature classes"
    outFlds = ["SHAPE@"] + [f.name for f in
                            arcpy.ListFields (outFc)
                            if f.type not in
                            ("Geometry", "Guid", "OID")]
    inputRows = [None] * len (outFlds)
    with arcpy.da.InsertCursor (outFc, outFlds) as iCurs:
        for i, fc in enumerate (mergeFcs, 1):
            print i, numFcs, "copying rows from", os.path.basename (fc)
            inFlds = ["SHAPE@"] + [f.name for f
                                   in arcpy.ListFields (fc)
                                   if f.name in outFlds]
            with arcpy.da.SearchCursor (fc, inFlds) as sCurs:
                for inRow in sCurs:
                    row = list (inputRows)
                    for n, fld in enumerate (inFlds):
                        if inRow [n]:
                            if fld in stringFlds:
                                try: row [outFlds.index
                                          (fld)] = str (inRow [n])
                                except: row [outFlds.index
                                             (fld)] = inRow [n]
                                
                            else:
                                row [outFlds.index
                                     (fld)] = inRow [n]
                    row [0] = inRow [0]
                    iCurs.insertRow (row)
    return outFc
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Tuesday, February 23, 2016

ArcGIS Erase without an ArcGIS advanced license

After creating code for converting Polygons to Lines and Vertices to Points, I decided to try my hand at creating an Erase function that didn't require ArcGIS's advance license. The code I came up with can erase polygons from lines, lines from lines, and polygons from polygons. The inputs are:


  • inFc: The input feature class
  • eraseFc: The feature class whose features will be erased from the input feature class
  • outLoc: The output file geodatabase
  • outName: The final name of the output feature class
  • workspace: A file geodatabase workspace
  • notices: Set to True for printed progress information


The code:

def Erase (inFc, eraseFc, outLoc = "", outName = "", workspace = "in_memory", notices = False):
    import os
    addToToc = arcpy.env.addOutputsToMap
    arcpy.env.addOutputsToMap = False
    if not outName:
        outName = os.path.basename (inFc) + "_Erase"
    if not outLoc:
        catPath = arcpy.Describe (inFc).catalogPath
        outLoc = os.path.dirname (catPath)
    outFc = UniqueFileName (outLoc, outName)
    outLoc, outName = os.path.split (outFc)

    garbage = []

    inType = arcpy.Describe (inFc).shapeType
    eraseType = arcpy.Describe (eraseFc).shapeType

    #lines and lines
    if inType == "Polyline" and eraseType == "Polyline":
        #buffer both lines: flat, left
        inBuffFc = UniqueFileName (workspace)
        arcpy.Buffer_analysis (inFc, inBuffFc, ".01 FEET")
        garbage += [inBuffFc]
        eraseBuffFc = UniqueFileName (workspace)
        arcpy.Buffer_analysis (eraseFc, eraseBuffFc, ".01 FEET")
        garbage += [eraseBuffFc]
        
        #union input buffer with erase buffer
        unionFc = UniqueFileName (workspace)
        arcpy.Union_analysis ([inBuffFc, eraseBuffFc], unionFc)
        garbage += [unionFc]
        #select union output by erase buffer
        unionLyr = UniqueLayerName (unionFc)
        garbage += [unionLyr]
        arcpy.SelectLayerByLocation_management (unionLyr, "WITHIN", eraseBuffFc)
        #switch union selection
        arcpy.SelectLayerByAttribute_management (unionLyr, "SWITCH_SELECTION")
        #check for selection
        if not arcpy.Describe (unionLyr).FIDSet:
            arcpy.CreateFeatureclass_management (outLoc, outName, "POLYLINE", inFc, "", "", inFc)
            for trash in garbage:
                arcpy.Delete_management (trash)
            return outFc
        #clip line
        arcpy.env.addOutputsToMap = addToToc
        arcpy.Clip_analysis (inFc, unionLyr, outFc)
        

    #lines and polygons
    elif inType == "Polyline" and eraseType == "Polygon":
        #buffer both lines: flat, left
        inBuffFc = UniqueFileName (workspace)
        arcpy.Buffer_analysis (inFc, inBuffFc, ".01 FEET")
        garbage += [inBuffFc]
        
        #union input buffer with erase buffer
        unionFc = UniqueFileName (workspace)
        arcpy.Union_analysis ([inBuffFc, eraseFc], unionFc)
        garbage += [unionFc]
        #select union output by erase buffer
        unionLyr = UniqueLayerName (unionFc)
        garbage += [unionLyr]
        arcpy.SelectLayerByLocation_management (unionLyr, "WITHIN", eraseFc)
        #switch union selection
        arcpy.SelectLayerByAttribute_management (unionLyr, "SWITCH_SELECTION")
        #check for selection
        if not arcpy.Describe (unionLyr).FIDSet:
            arcpy.CreateFeatureclass_management (outLoc, outName, "POLYLINE", inFc, "", "", inFc)
            for trash in garbage:
                arcpy.Delete_management (trash)
            return outFc
        #clip line
        arcpy.env.addOutputsToMap = addToToc
        arcpy.Clip_analysis (inFc, unionLyr, outFc)

    #polygons and polygons
    elif inType == "Polygon" and eraseType == "Polygon":
        unionFc = UniqueFileName (outLoc, outName)
        arcpy.env.addOutputsToMap = addToToc
        arcpy.Union_analysis ([inFc, eraseFc], unionFc, "NO_FID")
        arcpy.env.addOutputsToMap = False
        unionLyr = UniqueLayerName (unionFc)
        garbage += [unionLyr]
        arcpy.SelectLayerByLocation_management (unionLyr, "WITHIN", eraseFc)
        arcpy.DeleteFeatures_management (unionLyr)

        inFlds = [f.name for f in arcpy.ListFields (inFc)]
        delFlds = [f.name for f in arcpy.ListFields (unionLyr) if not f.name in inFlds and not
                   f.type in ("Geometry", "OID") and f.required == False]
        arcpy.DeleteField_management (unionLyr, delFlds)
        arcpy.env.addOutputsToMap = addToToc

    for trash in garbage:
        arcpy.Delete_management (trash)    
    
    return outFc

def UniqueFileName(Location = "in_memory", Name = "file", Extension = ""):
    """
    returns an unused file name
    'Location' will be the file folder
    'Name' will be the file name plus numeric extension if needed
    'Extension' is the file extension

    similar to arcpy.CreateUniqueName
    """
    import arcpy
    import os

    if Extension:

        outName = os.path.join (Location, Name + "." + Extension)
    else:
        outName = os.path.join (Location, Name)
    i = 0
    while arcpy.Exists (outName):
        i += 1
        if Extension:
            outName = os.path.join (Location, Name + "_" + str(i) + "." + Extension)
        else:
            outName = os.path.join (Location, Name + "_" + str(i))
        
    return outName

def UniqueLayerName(inFC = None, sql = None):
    """
    returns an unused feature layer name
    """
    import arcpy
    import os
    lyrName = os.path.join ("lyr0")
    i = 0
    while arcpy.Exists (lyrName):
        i += 1
        lyrName = "lyr" + str(i)
    if inFC:
        arcpy.MakeFeatureLayer_management (inFC, lyrName, sql) 
    return lyrName
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Wednesday, February 3, 2016

ArcGIS polygon to lines and feature vertices to points without an ArcGIS advanced license

I've started taking a little time to create python functions, using standard license code, that replicate some of ArcGIS's tools that are only available with an advanced license. So far, I've created functions for the Polygon To Line and the Feature Vertices To Points tools. Each make use of cursors and geometries to get the job done. The two functions take an input feature class's full catalog path, a workspace (file geodatabase or in_memory), and the file name for the new feature class, and return the full path to the function's output, which has been generated in the indicated workspace. Both also make use of my UniqueFileName function mentioned in a previous blog.

Which tools should I code next?

The code:

#polygons to lines
def PolygonToLines (inFc, workspace, fcName):

    outFc = UniqueFileName (workspace, fcName)
    try:
        #let esri do it (advanced license)
        arcpy.FeatureToLine_management (inFc, outFc)
        return outFc
    except:
        pass
    #no advanced license
    array = arcpy.Array ()
    sr = arcpy.Describe (inFc).spatialReference
    outPath, outName = os.path.split (outFc)
    arcpy.CreateFeatureclass_management (outPath,
                                         outName,
                                         "POLYLINE",
                                         spatial_reference = sr)
    with arcpy.da.InsertCursor (outFc, "SHAPE@") as iCurs:
        with arcpy.da.SearchCursor (inFc, "SHAPE@") as sCurs:
            for geom, in sCurs:
                for i in range (geom.partCount):
                    lastPart = None
                    part = geom.getPart (i)
                    for pnt in part:
                        if not pnt:
                            lastPart = None
                            continue
                        x = pnt.X
                        y = pnt.Y
                        if not lastPart:
                            lastPart = (x, y)
                            continue
                        thisPart = (x, y)
                        lastPoint = arcpy.Point ()
                        lastPoint.X = lastPart [0]
                        lastPoint.Y = lastPart [1]
                        thisPoint = arcpy.Point ()
                        thisPoint.X = thisPart [0]
                        thisPoint.Y = thisPart [1]
                        array.add (lastPoint)
                        array.add (thisPoint)
                        polyline = arcpy.Polyline (array)
                        array.removeAll ()
                        lastPart = thisPart
                        row = (polyline,)
                        iCurs.insertRow (row)
    del iCurs
    del sCurs
    return outFc


#feature vertices to points
def VerticesToPoints (inFc, workspace, fcName):

    outFc = UniqueFileName(workspace, fcName)
    try:
        #let esri do it (advanced license)
        arcpy.FeatureVerticesToPoints_management (inFc, outFc)
        return outFc
    except:
        pass
    #no advanced license
    sr = arcpy.Describe (inFc).spatialReference
    outPath, outName = os.path.split (outFc)
    arcpy.CreateFeatureclass_management (outPath,
                                         outName,
                                         "POINT",
                                         spatial_reference = sr)
    with arcpy.da.InsertCursor (outFc, "SHAPE@") as iCurs:
        with arcpy.da.SearchCursor (inFc, "SHAPE@") as sCurs:
            for geom, in sCurs:
                for i in range (geom.partCount):
                    part = geom.getPart (i)
                    for pnt in part:
                        if not pnt:
                            continue
                        row = (pnt,)
                        iCurs.insertRow (row)
    del iCurs
    del sCurs
    return outFc


#generate a unique ArcGIS file name
def UniqueFileName(location = "in_memory",
                   name = "file",
                   extension = ""):
    
    if extension:
        outName = os.path.join (location, name + "." + extension)
    else:
        outName = os.path.join (location, name)
    i = 0
    while arcpy.Exists (outName):
        i += 1
        if extension:
            outName = os.path.join (location,
                                    "{0}_{1}.{2}".format (name, i,
                                                          extension))
        else:
            outName = os.path.join (location,
                                    "{0}_{1}".format (name, i))
    return outName
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Thursday, December 17, 2015

Creating Unique Feature Class, Feature Layer, and Field Names with ArcPy

Most Python/ArcPy scripts frequently require the creation of new feature classes, feature layers, and, at times, fields. I developed scripts that return unique names for these so that I don't have to worry about my script throwing an error because files already exist, or worry that my script copies over existing data.

Unique file or feature class name:
def UniqueFileName(Location = "in_memory", Name = "file", Extension = ""):
    """
    returns an unused file name
    'Location' will be the file folder
    'Name' will be the file name plus numeric extension if needed
    'Extension' is the file extension

    similar to arcpy.CreateUniqueName
    """
    import arcpy
    import os

    if Extension:

        outName = os.path.join (Location, Name + "." + Extension)
    else:
        outName = os.path.join (Location, Name)
    i = 0
    while arcpy.Exists (outName):
        i += 1
        if Extension:
            outName = os.path.join (Location, 
                                    Name + "_" + str(i) + "." + Extension)
        else:
            outName = os.path.join (Location, Name + "_" + str(i))
        
    return outName
Unique feature layer name:
def UniqueLayerName(inFC = None, sql = None):
    """
    returns an unused feature layer name
    """
    import arcpy
    lyrName = "lyr0"
    i = 0
    while arcpy.Exists (lyrName):
        i += 1
        lyrName = "lyr" + str(i)
    if inFC:
        arcpy.MakeFeatureLayer_management (inFC, lyrName, sql) 
    return lyrName
Unique field name:
def UniqueFieldName(inFC, inField = "NEWFIELD"):
    """
    returns an unused field name
    """
    import arcpy

    Type = arcpy.Describe (inFC).dataType

    if Type == "ShapeFile":
        inField = inField[:10]

    fields = [f.name for f in arcpy.ListFields (inFC)]
    i = 0
    outField = inField
    while outField in fields:
        i += 1
        outField = inField + "_" + str(i)

        q = 0
        while Type == "ShapeFile" and len (outField) > 10:
            q -= 1
            outField = inField[:q] + "_" + str(i)

    return outField
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals