Saturday, October 10, 2015

Working Around Oracle's SQL Where Clause 1000 Item Limit with ArcPy

One option for working in an ArcSDE environment is to utilize Oracle databases. Oracle databases has an unfortunate limitation: they only allow for up to 1000 items when an SQL where clause is applied to them. When I create ArcGIS Python scripts, I will oftentimes perform a series of geoprocesses and throw field values (GUIDs perhaps) into a python list for features that meet certain criteria. I'll then want to filter my analyzed feature class by these values by means of a selection or through the creation of a feature layer. Both ArcPy's SelectLayerByAttribute_management and MakeFeatureLayer_management allow for the inclusion of an SQL where clause. In a large feature class, however, the 1000 item limit can be easily breached. To work around this limitation, I created a function that creates a feature layer from a list of values, even if the number of values is greater than 1000. The code first creates a new feature layer, and then performs multiple selections in 1000 item chunks. Once all the desired features have been selected, a final output feature layer is made and returned based on the selected features.

The code:

import arcpy

def LayerFromList (inLyr, inField, inList):

    """
    Returns an arcpy feature layer of the 
     features matching the values in 
     the inList in field inField
    inLyr can be a feature layer or 
     feature class
    inField is the field name in the 
     inLyr used for selection
    inList is a python list of values 
     to match to values in inField of inLyr
    Created to work around 1000 item 
     limit in Oracle database SQL querry
    """

    #convert list to strings
    mapList = map (str, inList)

    #Get Field Type
    fldType = [field.type for field
               in arcpy.ListFields (inLyr)
               if field.name == inField][0]
    
    #Make temporary feature layer to apply selections to
    tempLyr = "lyr0"
    num = 0
    #Get available layer name
    while arcpy.Exists (tempLyr):
        num += 1
        tempLyr = "lyr" + str (num)
    #Make layer
    arcpy.MakeFeatureLayer_management (inLyr, tempLyr)

    #features selected counter
    selected = 0
    #iterate
    while True:
        #Get next 1000 rows
        partial = mapList[selected:selected + 1000]
        #Exit while loop if there are no more items in the list
        if not partial:
            break

        #Create the SQL where cluase
        if fldType in ["String", "Guid"]:
            joinStr = "'{0}'".format ("', '".join (partial))
        else:
            joinStr = ", ".join (partial)
        sql = "{0} IN ({1})".format (arcpy.AddFieldDelimiters
                                     (tempLyr, inField), joinStr)

        #Add next 1000 rows to selection
        arcpy.SelectLayerByAttribute_management (tempLyr, 
                                                 "ADD_TO_SELECTION", 
                                                 sql)
        #Add 1000 to selected varaible
        selected += 1000


    #Make output feature layer
    num += 1
    outLyr = "lyr" + str (num)
    while arcpy.Exists (outLyr):
        num += 1
        outLyr = "lyr" + str (num)
    arcpy.MakeFeatureLayer_management (tempLyr, outLyr)

    #Delete temporary layer
    arcpy.Delete_management (tempLyr)
    
    return outLyr
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Tuesday, October 6, 2015

Zipping Shapefiles with Python

My previous employer liked using GIS Pro to collect GIS field data on iPads. My employer wanted to equip surveyors with a variety of GIS information to be used in the field - parcels, political boundaries, company assets, and others. To do so, we first had to create shapefiles and add them to zip files for our data to be compatible with GIS Pro. All this information required geoprocessing to make sure the right data was added to each surveyor's iPad without bogging down GIS Pro's software with a crippling overabundance of data. Once that processing was complete I needed to create a function that would get my shapefiles zipped up nicely and ready for upload. To start, I hunted down all the possible shapefile file extensions. With use of Python's zipfile module I developed my function to add all shapefile files to a new zipfile that inherits its file name from the shapefile. I also added the option to keep or delete the shapefile from its directory after zipping by setting the Boolean Delete variable to either True or False.

The code:
import os
import zipfile
import arcpy

def ZipShp (inShp, Delete = True):

    """
    Creates a zip file containing the input shapefile
    inputs -
    inShp: Full path to shapefile to be zipped
    Delete: Set to True to delete shapefile files after zip
    """
    
    #List of shapefile file extensions
    extensions = [".shp",
                  ".shx",
                  ".dbf",
                  ".sbn",
                  ".sbx",
                  ".fbn",
                  ".fbx",
                  ".ain",
                  ".aih",
                  ".atx",
                  ".ixs",
                  ".mxs",
                  ".prj",
                  ".xml",
                  ".cpg"]

    #Directory of shapefile
    inLocation = arcpy.Describe (inShp).path
    #Base name of shapefile
    inName = arcpy.Describe (inShp).baseName
    #Create zipfile name
    zipfl = os.path.join (inLocation, inName + ".zip")
    #Create zipfile object
    ZIP = zipfile.ZipFile (zipfl, "w")
    
    #Iterate files in shapefile directory
    for fl in os.listdir (inLocation):
        #Iterate extensions
        for extension in extensions:
            #Check if file is shapefile file
            if fl == inName + extension:
                #Get full path of file
                inFile = os.path.join (inLocation, fl)
                #Add file to zipfile
                ZIP.write (inFile, fl)
                break

    #Delete shapefile if indicated
    if Delete == True:
        arcpy.Delete_management (inShp)

    #Close zipfile object
    ZIP.close()

    #Return zipfile full path
    return zipfl

And for those of you without the ArcPy module:

import os
import zipfile

def ZipShp (inShp, Delete = True):

    """
    Creates a zip file containing the input shapefile
    inputs -
    inShp: Full path to shapefile to be zipped
    Delete: Set to True to delete shapefile files after zip
    """
    
    #List of shapefile file extensions
    extensions = [".shp",
                  ".shx",
                  ".dbf",
                  ".sbn",
                  ".sbx",
                  ".fbn",
                  ".fbx",
                  ".ain",
                  ".aih",
                  ".atx",
                  ".ixs",
                  ".mxs",
                  ".prj",
                  ".xml",
                  ".cpg",
                  ".shp.xml"]

    #Directory of shapefile
    inLocation = os.path.dirname (inShp)
    #Base name of shapefile
    inName = os.path.basename (os.path.splitext (inShp)[0])
    #Create zipfile name
    zipfl = os.path.join (inLocation, inName + ".zip")
    #Create zipfile object
    ZIP = zipfile.ZipFile (zipfl, "w")

    #Empty list to store files to delete
    delFiles = []
    
    #Iterate files in shapefile directory
    for fl in os.listdir (inLocation):
        #Iterate extensions
        for extension in extensions:
            #Check if file is shapefile file
            if fl == inName + extension:
                #Get full path of file
                inFile = os.path.join (inLocation, fl)
                #Add file to delete files list
                delFiles += [inFile]
                #Add file to zipfile
                ZIP.write (inFile, fl)
                break

    #Delete shapefile if indicated
    if Delete == True:
        for fl in delFiles:
            os.remove (fl)

    #Close zipfile object
    ZIP.close()

    #Return zipfile full path
    return zipfl
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals

Sunday, September 27, 2015

Determining Lengths of Polyline Features Regardless of Spatial Reference using ArcPy

Edit 2017: After all the below work, I've come to learn that there's a much easier way to determining feature lengths regardless of spatial reference with use of the field calcutor. This blog explains how. 

Here's an example:

arcpy.CalculateField_management ("InputLayer", "LenMiles", "!shape.length@miles!", "PYTHON_9.3")

Original Blog:

When working with GIS, knowing lengths of polyline features can be fundamental to accomplishing your goals. However, something so basic isn't necessarily always simple. I wanted to create a python function that would take any feature class and return each feature's length in a standard unit (meters). I wanted the code to work for any feature class, regardless of projection or how large or small an area of the globe the feature class's features covered.

Determining lengths of features in a projected feature class is easy. ArcPy allows you to do so by the feature class's spatial reference object. This object has the metersPerUnit method, which returns the conversion factor for the feature class's units to meters. Once this is known, a cursor can be used to obtain the length of each feature. This distance can be multiplied by the conversion factor for the standardized distance.

Geographic coordinate systems do not have a measure unit, so determining each feature's length is not as straight forward a task. Initially I couldn't find a solution, and so asked this question on GIS StackExchange. The answer presented the solution of setting the feature class's coordinate system within a cursor. Picking which coordinate system also wasn't easy. The answer suggested WGS Web Mercator. However, one of the commenters warned that WGS Web Mercator isn't actually in meters. He suggested using the UTM zone the data fell in as the projection. This was a good suggestion, except for the fact that I wanted my code to be robust enough to work regardless of where it was geographically. I had to find a way to determine which UTM zone each feature fell within. With a little poking around I came across Calculate UTM Zone. This tool appears to be designed for use with data driven pages, but I found it suited my needs relatively well. The output, which populates a field with WKT, wasn't perfect however. I asked this question for suggestions on how I could determine each feature's UTM zone from its output consistently, but didn't find an answer. Initially I came up with my own string manipulation method, but later discovered the loadFromString () method of the spatial reference object, that lets you update a spatial reference from well known text. From here, it was simply a matter of projecting each feature into its corresponding UTM zone and creating the dictionary.

The code:

def FeatureLengths (inFc):
 
    """
    input - inFc: the input feature class
    returns dictionary -
    key: feature OID
    value: distance in meters
    """
 
    #get spatial reference type
    sr = arcpy.Describe (inFc).spatialReference
    srType = sr.type
 
 
    #Projected coordinate systems
    if srType == "Projected":
        #meters per unit
        metersPerUnit = sr.metersPerUnit
 
        #out Di
        return dict ([(oid, dist * metersPerUnit) for
                      oid, dist in
                      arcpy.da.SearchCursor
                      (inFc, ["OID@", "SHAPE@LENGTH"])])
     
    #Geographic coordinate systems
    #add UTM field to feature class
    flds = [f.name for f in arcpy.ListFields (inFc)]
    utmFld = "UTM"
    i = 0
    #Make sure field name is not in use
    while utmFld in flds:
        utmFld = "UTM{0}".format (i)
        i += 1
    arcpy.AddField_management (inFc,
                               utmFld,
                               "TEXT",
                               field_length = 1000)
 
    
    #calculate UTM zone
    arcpy.CalculateUTMZone_cartography (inFc, utmFld)

    #add sql field delimiters
    delimUtmFld = arcpy.AddFieldDelimiters (inFc, utmFld)

    #output dictionary that will hold oids and distances in meters
    outDi = {}
    #list to hold utm zones
    utms = []

    #get unique layer name
    layer = "lyr"
    while arcpy.Exists (layer):
        i += 1
        layer = "lyr" + str (i)
 
    #access table
    with arcpy.da.SearchCursor (inFc, [utmFld]) as cursor:
        for utmWkt, in cursor:
            #check if utm has been iterated
            if utmWkt in utms:
                continue
            utms += [utmWkt]

            #create new spatial reference
            sr = arcpy.SpatialReference ()

            #update spatial reference with wkt
            sr.loadFromString (utmWkt)

            #create sql to filter for utm zone
            sql = "{0} = '{1}'".format (delimUtmFld, utmWkt)

            #create feature layer per utm zone
            arcpy.MakeFeatureLayer_management (inFc, layer, sql)

            #iterate layer with projection
            with arcpy.da.SearchCursor (layer,
                                        ["OID@", "SHAPE@LENGTH"],
                                        spatial_reference
                                        = sr) as lyrCursor:
                for oid, shpLen in lyrCursor:
                    #update dictionary with length
                    if shpLen:
                        outDi [oid] = shpLen
                    #if shape error, return zero
                    else:
                        outDi [oid] = 0
            del lyrCursor

            #Delete layer
            arcpy.Delete_management (layer)
            
    del cursor

    #delete utm zone field
    arcpy.DeleteField_management (inFc, utmFld)
    
    #return out dictionary
    return outDi
profile for Emil Brundage at Geographic Information Systems Stack Exchange, Q&A for cartographers, geographers and GIS professionals