ArcToolbox Tool: Add [ST_CON_ABR] to Hennepin County, MN Centerlines

One of the great advancements over the last decade plus in GIS is that government agencies have started to move away from a “recover-our-cost” mentality to more of an “Open Data”. Minnesota, for example, has launched their Geospatial Commons as a platform for sharing data.

And while getting free, authoritative data is awesome, it can leave you in a bind if the structure of the data changes. Sometime between April and September, Hennepin County, Minnesota, changed the schema of their publicly available street centerlines data.

The data used to have both a full, concatenated street name field ([ST_CONCAT]) and an abbreviated version ([ST_CON_ABR]). A record might have “James Lofton Avenue North” and “James Lofton Ave N”, respectively, in these two fields. The abbreviated version was nice for labeling but it disappeared from the most recent updates.

So, as you  might guess from the fact that I’m posting about it, I wrote a script to add that field in and launch it from an ArcToolbox Tool. Nothing fancy going on in the code, just a series of replaces, depending on the field. Using dictionaries instead of arrays of paired values might have been better but the script takes just a few seconds to run so I can live with it as-is.

The list of street type abbreviations came from a combination of ESRI’s standards and those found in an older version of the Hennepin County data. There were no conflicting abbreviations between the two. The code warns if a street name occurs in the data that is not in the list.

While I’m including the code here for reference, it’s probably best to download the code from GitHub.

#-------------------------------------------------------------------------------
# Name:        usi_dataprep_Add_STCONABR
#
# Purpose:     This can be used to add [CT_CON_ABR] to Hennepin County, MN
#              centerlines. This is a concatenated, abbreviated full name of
#              the street. This used to be included in the data but
#              disappeared from the downloads in the summer of 2017.
#
#              Data available at: http://www.hennepin.us/gisopendata
#
# Author:      mrantala
#
# Created:     2017.10.04
#
#------------------------------------------------------------------------------

import arcpy

############################################
## Custom Variables
#These are the fields that are concatenated. Hennepin has others but were always blank.
requiredFieldList = ["ST_PRE_DIR","ST_PRE_TYP","ST_NAME","ST_POS_TYP","ST_POS_DIR"]
#This is the name of the field to add
newFieldName = "ST_CON_ABR"

#These are the abbreviations for [ST_POS_TYPE]. The list was created using a sample of
#Hennepin's centerline data & Esri Tech article: http://support.esri.com/en/technical-article/000008454
# Note that I intentionally left cases in where there is no abbreviation (Fall, for example) as a means of
#documenting the fact that it should NOT change.

abbList = []
abbList.append(["Alcove","Alcove"]) #Hennepin Specific
abbList.append(["Alley","Aly"])
abbList.append(["Annex","Anx"])
abbList.append(["Arcade","Arc"])
abbList.append(["Avenue","Ave"])
abbList.append(["Bay","Bay"]) #Hennepin Specific
abbList.append(["Bayoo","Byu"])
abbList.append(["Beach","Bch"])
abbList.append(["Bend","Bnd"])
abbList.append(["Bluff","Blf"])
abbList.append(["Bluffs","Blfs"])
abbList.append(["Bottom","Btm"])
abbList.append(["Boulevard","Blvd"])
abbList.append(["Branch","Br"])
abbList.append(["Bridge","Brg"])
abbList.append(["Brook","Brk"])
abbList.append(["Brooks","Brks"])
abbList.append(["Burg","Bg"])
abbList.append(["Burgs","Bgs"])
abbList.append(["Bypass","Byp"])
abbList.append(["Camp","Cp"])
abbList.append(["Canyon","Cyn"])
abbList.append(["Cape","Cpe"])
abbList.append(["Causeway","Cswy"])
abbList.append(["Center","Ctr"])
abbList.append(["Centers","Ctrs"])
abbList.append(["Crossings","Crossings"]) #Hennepin Specific
abbList.append(["Crossroad","Xrd"])
abbList.append(["Chase","Chase"]) #Hennepin Specific
abbList.append(["Circle","Cir"])
abbList.append(["Circles","Cirs"])
abbList.append(["Cliff","Clf"])
abbList.append(["Cliffs","Clfs"])
abbList.append(["Club","Clb"])
abbList.append(["Close","Close"]) #Hennepin Specific
abbList.append(["Common","Cmn"])
abbList.append(["Commons","Cmns"]) #Hennepin Specific
abbList.append(["Corner","Cor"])
abbList.append(["Corners","Cors"])
abbList.append(["Corridor","Corridor"]) #Hennepin Specific
abbList.append(["Course","Crse"])
abbList.append(["Court","Ct"])
abbList.append(["Courts","Cts"])
abbList.append(["Cove","Cv"])
abbList.append(["Coves","Cvs"])
abbList.append(["Creek","Crk"])
abbList.append(["Crescent","Cres"])
abbList.append(["Crest","Crst"])
abbList.append(["Cross","Cross"]) #Hennepin Specific
abbList.append(["Crossing","Xing"])
abbList.append(["Curve","Curve"])
abbList.append(["Dale","Dl"])
abbList.append(["Dam","Dm"])
abbList.append(["Divide","Dv"])
abbList.append(["Down","Down"]) #Hennepin Specific
abbList.append(["Downs","Downs"]) #Hennepin Specific
abbList.append(["Drive","Dr"])
abbList.append(["Drives","Drs"])
abbList.append(["Edge","Edge"]) #Hennepin Specific
abbList.append(["Entry","Entry"]) #Hennepin Specific
abbList.append(["Estate","Est"])
abbList.append(["Estates","Ests"])
abbList.append(["Expressway","Expy"])
abbList.append(["Extension","Ext"])
abbList.append(["Extensions","Exts"])
abbList.append(["Fall","Fall"])
abbList.append(["Falls","Fls"])
abbList.append(["Ferry","Fry"])
abbList.append(["Field","Fld"])
abbList.append(["Fields","Flds"])
abbList.append(["Flat","Flt"])
abbList.append(["Flats","Flts"])
abbList.append(["Ford","Frd"])
abbList.append(["Fords","Frds"])
abbList.append(["Forest","Frst"])
abbList.append(["Forge","Frg"])
abbList.append(["Forges","Frgs"])
abbList.append(["Fork","Frk"])
abbList.append(["Forks","Frks"])
abbList.append(["Fort","Ft"])
abbList.append(["Freeway","Fwy"])
abbList.append(["Gables","Gables"]) #Hennepin Specific
abbList.append(["Garden","Gdn"])
abbList.append(["Gardens","Gdns"])
abbList.append(["Gate","Gate"]) #Hennepin Specific
abbList.append(["Gateway","Gtwy"])
abbList.append(["Glade","Glade"]) #Hennepin Specific
abbList.append(["Glen","Gln"])
abbList.append(["Glens","Glns"])
abbList.append(["Green","Grn"])
abbList.append(["Greens","Grns"])
abbList.append(["Greenway","Greenway"]) #Hennepin Specific
abbList.append(["Grove","Grv"])
abbList.append(["Groves","Grvs"])
abbList.append(["Harbor","Hbr"])
abbList.append(["Harbors","Hbrs"])
abbList.append(["Haven","Hvn"])
abbList.append(["Heights","Hts"])
abbList.append(["Highway","Hwy"])
abbList.append(["Hill","Hl"])
abbList.append(["Hills","Hls"])
abbList.append(["Hollow","Holw"])
abbList.append(["Horn","Horn"]) #Hennepin Specific
abbList.append(["Inlet","Inlt"])
abbList.append(["Island","Is"])
abbList.append(["Islands","Iss"])
abbList.append(["Isle","Isle"])
abbList.append(["Junction","Jct"])
abbList.append(["Junctions","Jcts"])
abbList.append(["Key","Ky"])
abbList.append(["Keys","Kys"])
abbList.append(["Knoll","Knl"])
abbList.append(["Knolls","Knls"])
abbList.append(["Lake","Lk"])
abbList.append(["Lakes","Lks"])
abbList.append(["Land","Land"])
abbList.append(["Landing","Lndg"])
abbList.append(["Lane","Ln"])
abbList.append(["Light","Lgt"])
abbList.append(["Lights","Lgts"])
abbList.append(["Loaf","Lf"])
abbList.append(["Lock","Lck"])
abbList.append(["Locks","Lcks"])
abbList.append(["Lodge","Ldg"])
abbList.append(["Loop","Loop"])
abbList.append(["Mall","Mall"])
abbList.append(["Manor","Mnr"])
abbList.append(["Manors","Mnrs"])
abbList.append(["Meadow","Mdw"])
abbList.append(["Meadows","Mdws"])
abbList.append(["Mews","Mews"])
abbList.append(["Mill","Ml"])
abbList.append(["Mills","Mls"])
abbList.append(["Mission","Msn"])
abbList.append(["Motorway","Mtwy"])
abbList.append(["Mount","Mt"])
abbList.append(["Mountain","Mtn"])
abbList.append(["Mountains","Mtns"])
abbList.append(["Neck","Nck"])
abbList.append(["Orchard","Orch"])
abbList.append(["Oval","Oval"])
abbList.append(["Overpass","Opas"])
abbList.append(["Park","Park"])
abbList.append(["Parks","Park"])
abbList.append(["Parkway","Pkwy"])
abbList.append(["Parkways","Pkwy"])
abbList.append(["Pass","Pass"])
abbList.append(["Passage","Psge"])
abbList.append(["Path","Path"])
abbList.append(["Pike","Pike"])
abbList.append(["Pine","Pne"])
abbList.append(["Pines","Pnes"])
abbList.append(["Place","Pl"])
abbList.append(["Plain","Pln"])
abbList.append(["Plains","Plns"])
abbList.append(["Plaza","Plz"])
abbList.append(["Point","Pt"])
abbList.append(["Points","Pts"])
abbList.append(["Port","Prt"])
abbList.append(["Ports","Prts"])
abbList.append(["Prairie","Pr"])
abbList.append(["Radial","Radl"])
abbList.append(["Railroad","Railroad"]) #Hennepin Specific
abbList.append(["Ramp","Ramp"])
abbList.append(["Ranch","Rnch"])
abbList.append(["Rapid","Rpd"])
abbList.append(["Rapids","Rpds"])
abbList.append(["Rest","Rst"])
abbList.append(["Ridge","Rdg"])
abbList.append(["Ridges","Rdgs"])
abbList.append(["Rise","Rise"]) #Hennepin Specific
abbList.append(["River","Riv"])
abbList.append(["Road","Rd"])
abbList.append(["Roads","Rds"])
abbList.append(["Route","Rte"])
abbList.append(["Row","Row"])
abbList.append(["Rue","Rue"])
abbList.append(["Run","Run"])
abbList.append(["Shoal","Shl"])
abbList.append(["Shoals","Shls"])
abbList.append(["Shore","Shr"])
abbList.append(["Shores","Shrs"])
abbList.append(["Skies","Skies"]) #Hennepin Specific
abbList.append(["Skyway","Skwy"])
abbList.append(["Spring","Spg"])
abbList.append(["Springs","Spgs"])
abbList.append(["Spur","Spur"])
abbList.append(["Spurs","Spur"])
abbList.append(["Square","Sq"])
abbList.append(["Squares","Sqrs"])
abbList.append(["Station","Sta"])
abbList.append(["Stravenue","Stra"])
abbList.append(["Stream","Strm"])
abbList.append(["Street","St"])
abbList.append(["Streets","Sts"])
abbList.append(["Summit","Smt"])
abbList.append(["Terrace","Ter"])
abbList.append(["Throughway","Trwy"])
abbList.append(["Trace","Trce"])
abbList.append(["Track","Trak"])
abbList.append(["Trafficway","Trfy"])
abbList.append(["Trail","Trl"])
abbList.append(["Tunnel","Tunl"])
abbList.append(["Turn","Turn"]) #Hennepin Specific
abbList.append(["Turnpike","Tpke"])
abbList.append(["Underpass","Upas"])
abbList.append(["Union","Un"])
abbList.append(["Unions","Uns"])
abbList.append(["Valley","Vly"])
abbList.append(["Valleys","Vlys"])
abbList.append(["Viaduct","Via"])
abbList.append(["View","Vw"])
abbList.append(["Views","Vws"])
abbList.append(["Village","Vlg"])
abbList.append(["Villages","Vlgs"])
abbList.append(["Ville","Vl"])
abbList.append(["Vista","Vis"])
abbList.append(["Walk","Walk"])
abbList.append(["Walks","Walk"])
abbList.append(["Wall","Wall"])
abbList.append(["Way","Way"])
abbList.append(["Ways","Ways"])
abbList.append(["Well","Wl"])
abbList.append(["Wells","Wls"])

#List of changes for [St_POS_Dir]
posDirList = [["North","N"],["East","E"],["South","S"],["West","W"],["Northeast","NE"],["Northwest","NW"],["Southeast","SE"],["Southwest","SW"]]
preDirList = [["North","N"],["East","E"],["South","S"],["West","W"]]
############################################
## Read Arguments

if (len(sys.argv) > 1):
    inFC = sys.argv[1]

############################################
# General Purpose Functions
def printit(inputString):
    try:
        print(inputString)
        arcpy.AddMessage(str(inputString))
    except:
        pass

def printerror(inputString):
    print (inputString)
    arcpy.AddError(inputString)

def getField(inFeatureClass, inFieldName):
  fieldList = arcpy.ListFields(inFeatureClass)
  for iField in fieldList:
    if iField.name.lower() == inFieldName.lower():
      return iField
  return None

def fieldExists(inFeatureClass, inFieldName):
  return getField(inFeatureClass,inFieldName) <> None

############################################
# Initial QC

def initialQC():
    if (arcpy.Exists(inFC)):
        printit("PASS: Feature Class {} Exists".format(inFC))
    else:
        printerror("ERROR: Feature Class {} Does Not Exist, Cancelling...".format(inFC))
        return False

    for iFld in requiredFieldList:
        if (fieldExists(inFC,iFld)):
            printit("PASS: Feature Class {} Has Field [{}]".format(inFC,iFld))
        else:
            printerror("ERROR: Feature Class {} Does Not Have Field [{}], Cancelling...".format(inFC,iFld))
            return False

    if not (fieldExists(inFC,newFieldName)):
        printit("GOOD: Feature Class {} Does Not Already Have Field [{}]".format(inFC,newFieldName))
        printit(" ADDING Field [{}]".format(newFieldName))
        try:
            arcpy.AddField_management(in_table=inFC, field_name=newFieldName, field_type="TEXT", field_precision="", field_scale="", field_length="100", field_alias="", field_is_nullable="NULLABLE", field_is_required="NON_REQUIRED", field_domain="")
        except:
            printerror("ERROR: Error While Adding Field [{}], Cancelling...".format(newFieldName))
            return False
        if not (fieldExists(inFC,newFieldName)):
            printerror("ERROR: Unable to Add Field [{}], Cancelling...".format(newFieldName))
            return False
    else:
        printerror("ERROR: Feature Class {} Already Has Field [{}], Cancelling...".format(inFC,newFieldName))
        return False

    return True

############################################
# Main

def makeSubstitution(inList,inValue,inFieldName):
    for iAbbreviationPr in inList:
        if (inValue == iAbbreviationPr[0]): #Found a Match
            return iAbbreviationPr[1]
    printit("WARNING: [{}] of {} does not have a value in the abbreviation list! Potential Error...".format(inFieldName,inValue))
    return inValue

def main():
    cursorFieldList = requiredFieldList
    cursorFieldList.append(newFieldName)

    try:
        iUCursor = arcpy.da.UpdateCursor(inFC,cursorFieldList)
        iRowCount = 0
        iRowMax = 1
        for uRow in iUCursor:

            #Just to give user an indicator that progress is being made
            if (iRowCount>iRowMax):
                printit(" {}".format(iRowCount))
                iRowMax *= 10
                iRowCount+=1

            abbreviateConcatenatedName = ""
            iFldIndex = 0
            for iFld in requiredFieldList:


                if (iFld == newFieldName):
                    uRow[iFldIndex] = abbreviateConcatenatedName
                    iUCursor.updateRow(uRow)
                else:
                    iValue = uRow[iFldIndex].strip() #Strip is just a safe-guard


                    if ((iValue != "") and (iValue != None)):
                        if (iFld == "ST_PRE_DIR"):
                            iValue= makeSubstitution(preDirList,iValue,"ST_PRE_DIR")
                        if (iFld == "ST_POS_TYP"):
                            iValue= makeSubstitution(abbList,iValue,"ST_POS_TYPE")
                        if (iFld == "ST_POS_DIR"):
                            iValue = makeSubstitution(posDirList,iValue,"ST_POS_DIR")

                        if (abbreviateConcatenatedName == ""):
                            abbreviateConcatenatedName = iValue
                        else:
                            abbreviateConcatenatedName+=" "+iValue
                iFldIndex += 1

        del iUCursor
    except RuntimeError as e:
        printerror("ERROR: Error {} Occurred, Cancelling...".format(e))
        try:
            del iUCursor
            del uRow
        except:
            return False
    return True


if __name__ == '__main__':
    if (initialQC() == True):
        if (main() == True):
            printit("Done!")

Converting MXD to Layer file in Arcpy

Working on doing some advanced ArcGIS server printing and had the need to batch convert many existing .mxd files to .lyr files. So instead of opening up X number of map documents, thought I would do it via code. All of my .mxds in this case had just one data frame so the process was pretty simple–I add an empty group layer (Thanks Petr Krebs for the idea), copy all the existing layers into it, and save it out as a layer file.

I created an ArcGIS toolbox with two options–one to convert a single .mxd and one to batch convert an entire folder. To use it, make sure to have the EmptyGroup.lyr in the same directory as the .py file.

Here is the raw code or git it:


import os
import arcpy
import inspect
import glob
import uuid
import inspect

codeDir = os.path.dirname(inspect.getfile(inspect.currentframe()))
EmptyGroupLayerFile = codeDir+"/EmptyGroup.lyr"
inArg1 = sys.argv[1]
inArg2 = sys.argv[2]

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

def makeLyrFromMXD(inMXD, outLyr):
    if not (os.path.exists(inMXD)):
        printit( "ERROR: {} does not exist".format(inMXD))
        return False
    if not (os.path.exists(EmptyGroupLayerFile)):
        printit( "ERROR: {} does not exist".format(EmptyGroupLayerFile))
        return False
    if  (os.path.exists(outLyr)):
        printit( "Skipping: {} already exists".format(outLyr))
        return True

    printit( "Making Layer file: {0}".format(outLyr))

    mxd = arcpy.mapping.MapDocument(inMXD)
    ###Right now, just doing the first Dataframe, this could be modified
    df = arcpy.mapping.ListDataFrames(mxd)[0]

    theUUID = str(uuid.uuid1())

    iGroupLayerRaw = arcpy.mapping.Layer(EmptyGroupLayerFile)
    iGroupLayerRaw.name = theUUID
    arcpy.mapping.AddLayer(df,iGroupLayerRaw,"TOP")
    groupBaseName = os.path.basename(outLyr).split(".")[0]

    for lyr in arcpy.mapping.ListLayers(df):
        if not (lyr.name == theUUID):
            if (lyr.longName == lyr.name):
                arcpy.mapping.AddLayerToGroup (df, iGroupLayer, lyr, "Bottom")
        else:
            iGroupLayer = lyr

    iGroupLayer.name = groupBaseName
    arcpy.SaveToLayerFile_management(iGroupLayer, outLyr)
    return os.path.exists(outLyr)

def doMultiple(inDir,outDir):
    for iMxd in glob.glob(inDir+"/*.mxd"):
        lyrFile = outDir+"/"+os.path.basename(iMxd).lower().replace(".mxd",".lyr")
        makeLyrFromMXD(iMxd, lyrFile)

if(not os.path.exists(EmptyGroupLayerFile)):
    printit("Error: {} is missing, can not run.".format(EmptyGroupLayerFile))
else:
    if (os.path.isdir(inArg1) and (os.path.isdir(inArg2))):
        doMultiple(inArg1,inArg2)
    elif (os.path.isfile(inArg1)):
        if (os.path.exists(inArg2)):
            printit("Error: {} already exists".format(inArg2))
        else:
            makeLyrFromMXD(inArg1,inArg2)
    else:
        printit("Unable to understand input parameters")

Calling os.startfile and webbrowser.open from ArcGIS.

Recently I’ve created Python add-ins for data entry for our staff. Most of these have a toolbar with a “Help” button that opens a help file in .pdf format.

Sample python add-in toolbar.
Sample python add-in toolbar.

The first add-in was for ArcCatalog and this worked splendidly. I was using os.startfile(path to help.pdf).

However, when I started doing ArcMap add-ins, clicking the Help button would open the help.pdf but ArcMap would crash. Oops!

Luckily the Python development team at Esri already had a blog post about this at their ArcPy Café blog.

They report that the root of the problem is “conflicts in the way the Windows libraries expect to be called, they can fail or crash when called within ArcGIS for Desktop in an add-in script or geoprocessing script tool”. But this can be overcome by using a decorator function that calls os.startfile from a new thread. Another function effected by these conflicts is webbrowser.open.

Example code is shown below:

import functools
import os
import threading
import webbrowser
 
# A decorator that will run its wrapped function in a new thread
def run_in_other_thread(function):
    # functool.wraps will copy over the docstring and some other metadata
    # from the original function
    @functools.wraps(function)
    def fn_(*args, **kwargs):
        thread = threading.Thread(target=function, args=args, kwargs=kwargs)
        thread.start()
        thread.join()
    return fn_
 
# Our new wrapped versions of os.startfile and webbrowser.open
startfile = run_in_other_thread(os.startfile)
openbrowser = run_in_other_thread(webbrowser.open)

Then whenever you call startfile or openbrowser, it will be routed through your decorator function and, as far as I’ve been able to tell, works fine without crashing your ArcMap session.

Cheers!

Quick & Dirty ArcPy: Listing Data Sources

I just had the need to go through a directory containing many (100+) layer files (.lyr) and verify the data sources in each. I could have loaded each into ArcMap and checked the properties, but choose not to. Here’s the bare-bones script I used instead:

import arcpy, glob,os

theDir = r"L:\gdrs\data\org\us_mn_state_dnr\elev_minnesota_lidar\\"
os.chdir(theDir)

for iFile in glob.glob("*.lyr"):
    print iFile
    lyr = arcpy.mapping.Layer(iFile)
    for i in arcpy.mapping.ListLayers(lyr):
        try:
            print "    {0}: {1}".format(i,i.dataSource)
        except:
            print "    {0}: Does not support dataSource".format(i)

print "Done!"

Garmin GPX to Shapefile (SHP) conversion GPX2Shp.py

I mentioned using Tapiriik to batch download my entire Garmin Connect history–over 1,000 separate .GPX files. I found several tools to convert .GPX to shapefiles that worked but none seemed to recognize my heart rate data.

The trick is Garmin extends the GPX specification to incorporate the heart rate:

xmlns:gpxtpx="http://www.garmin.com/xmlschemas/TrackPointExtension/v1"

Each track point looks like this:

     <trkpt lat=”43.68346489146352″ lon=”-92.99583793468773″>
        
        <ele>296.20001220703125
        <extensions>
          <gpxtpx:TrackPointExtension>
            <gpxtpx:hr>86
          gpxtpx:TrackPointExtension>
        </extensions>
      trkpt>

 

Since the first few exiting GPX converters failed to meet my needs, I decided to make my own, at least partially.

I used Joel Lawhead of GeospatialPython.com‘s pyshp library to handle writing the shapefile. I added some basic loop and I stuck a template.prj in the workspace that gets copied once for each shapefile.

Otherwise, nothing too fancy going on.  The code can be downloaded from Github.

import glob, os
import xml.etree.ElementTree as ET
import shapefile
import shutil

theNS = "{http://www.topografix.com/GPX/1/1}".lower()
theNS2 = "{http://www.garmin.com/xmlschemas/TrackPointExtension/v1}".lower()
templatePRJfile = "template.prj"

def elementIs(inElement,inTag):
    item1 = inTag.lower()
    item2 = elementName(inElement)
    return (inTag.lower() == elementName(inElement).lower())

def elementName(inElement):
    item1= inElement.tag.lower().replace(theNS,"").replace(theNS2,"")
    return item1

def convertTimeToSeconds(inTime):
    theSeconds = -1

    if (inTime.count(":")) == 2:
        try:
            inHour = inTime.split(":")[0]
            inMin = inTime.split(":")[1]
            inSec = inTime.split(":")[2]

            totalSec = float(inSec)
            totalSec += (float(inMin) * 60)
            totalSec += (float(inHour) * 3600)
            theSeconds = totalSec
        except:
            pass

    return theSeconds


def writeSHP(inSourceFile,inTrkList):
    w = shapefile.Writer(shapefile.POINT)
    w.field("file")
    w.field("segment","N","8",0)
    w.field("vertex","N","8",0)
    w.field("datetime","C",30)
    w.field("date","C","10",0)
    w.field("time","C","8",0)
    w.field("sec","N","8",0)
    w.field("isec","N","8",0)
    w.field("totsec","N","8",0)
    w.field("elev","N","24",14)
    w.field("hr","N","8",0)
    w.field("last","N","1",0)
    w.field("lat","N","24",16)
    w.field("lon","N","24",16)

    iTrkSegIndex = 0
    startSec =-1
    prevSec = -1
    for iTrkSeg in inTrkList:
        iTrkPtIndex = 0
        for iTrkPtDict in iTrkSeg:
            thisLine = "{0},{1},{2},*time*,*ele*,*hr*,*lat*,*lon*".format(inSourceFile,iTrkSegIndex,iTrkPtIndex)

            theLat = None
            if (iTrkPtDict.has_key('lat')):
                try:
                    theLat = float(iTrkPtDict['lat'])
                except:
                    pass

            theLon = None

            if (iTrkPtDict.has_key('lon')):
                try:
                    theLon = float(iTrkPtDict['lon'])
                except:
                    pass

            theDate = None
            theTime = None
            theSeconds = -1
            segSeconds = -1
            totSeconds = -1

            if (iTrkPtDict.has_key('time')):
                theDateTime = iTrkPtDict['time']
                if ("T" in theDateTime):
                    theDate = theDateTime.split("T")[0]
                    theTimePlue = theDateTime.split("T")[1]
                    if ("+" in theTimePlue):
                        theTime = theTimePlue.split("+")[0]
                        theSeconds = convertTimeToSeconds(theTime)

                        if (prevSec < 0):
                            prevSec = theSeconds
                        if (startSec<0):
                            startSec = theSeconds

                        segSeconds = theSeconds - prevSec
                        prevSec = theSeconds
                        totSeconds = theSeconds - startSec
            else:
                theDateTime = None

            if (iTrkPtDict.has_key('ele')):
                theElev = iTrkPtDict['ele']
            else:
                theElev = None

            if (iTrkPtDict.has_key('hr')):
                theHR = iTrkPtDict['hr']
            else:
                theHR = None

            if (iTrkPtIndex == len(iTrkSeg) - 1):
                theLast = 1
            else:
                theLast = 0

            w.point(theLon, theLat)
            try:
                                  w.record(inSourceFile.replace(".gpx",""),iTrkSegIndex,iTrkPtIndex,theDateTime,theDate,theTime,theSeconds,segSeconds,totSeconds,theElev,theHR,theLast,theLat,theLon)

            except:
                print "############## ERROR ####################"
            iTrkPtIndex+=1

        iTrkSegIndex+=1


    w.save(inSourceFile.lower().replace(".gpx",""))
    w = None
    if (os.path.exists(templatePRJfile)):
        newPRJFN = inSourceFile.lower().replace(".gpx",".prj")
        shutil.copyfile(templatePRJfile,newPRJFN)

def mainLoop():
    for iFile in glob.glob("*.gpx"):
        print iFile
        tree = ET.parse(iFile)
        root=tree.getroot()

        theTrkList = []

        for iRoot in root:
            if elementIs(iRoot,"trk"): #"http://www.topografix.com/gpx/1/1}trk" in iRoot.tag.lower():
                for iTrkSeg in iRoot:
                    if not elementIs(iTrkSeg,"trkseg"):
                        continue
                    thisTrk = []

                    pntIndex = 0
                    for iTrkPt in iTrkSeg:
                        if not elementIs(iTrkPt,"trkpt"):
                            continue
                        trkPntDict = dict()
                        trkPntDict["pntIndex"] = pntIndex
                        trkPntDict['lat'] = iTrkPt.get('lat')
                        trkPntDict['lon'] = iTrkPt.get('lon')

                        pntIndex+=1
                        for iElem in iTrkPt:
                            if elementIs(iElem,"extensions"):
                                for iSubElem in iElem:
                                    if (elementIs(iSubElem,"TrackPointExtension")):
                                        for iExtensionElem in iSubElem:
                                            if elementIs(iExtensionElem,"hr"):
                                                trkPntDict[elementName(iExtensionElem)] = iExtensionElem.text
                            else:
                                trkPntDict[elementName(iElem)] = iElem.text

                        #print trkPntDict
                        thisTrk.append(trkPntDict)

                    theTrkList.append(thisTrk)
        writeSHP(iFile.lower(), theTrkList)


theLineList = mainLoop()

 

Debugging a Python Scheduled Task

I have been working on a python script that I want (NEED) to run as a scheduled task on a remote machine.  I got to the point that the script did exactly what I needed when I was interactively running it in a Windows session but had problems when running it as a scheduled task.  The debugging process was cumbersome–make a change, schedule a task to run it, log out of the machine, and wait.  The log back in and repeat the process.

That got old.

So I wrote a script  (tester.py) that calls any other python scripts in the same directory that (1) start with “test_” and (2) there is not a corresponding file with the same base name and “.start” extension.  It would launch “test_BaBing.py” as long as there is not a “test_BaBing.start” in the same directory.  Tester.py continued to run, looping every 60 seconds, until tester.stop exists.

This made the process easier because I could work on my local machine, editing the problematic script, saving changes and within 60 seconds it would be launched on the remote machine.  I could view the results, make additional edits, delete the .start file and it would launch again within 60 seconds.

Within a couple minutes I was able to determine the problem (path related) and fix it.

Happy programmer.

<disclaimer>I would recommend using this only while debugging a script–routinely running it could be a security risk since someone could copy a destructive python script into the directory and this would run it.</disclaimer>

Download: tester.py

import sys, string, os
import glob
import datetime, shutil
import time, inspect
import getpass

totalstarttime = datetime.datetime.now()

dateString = datetime.date.today().strftime("%Y%m%d_")+datetime.datetime.now().strftime("%H%M%S") #datetime.date.today().strftime("%Y%m%d")
debugfile = inspect.getfile(inspect.currentframe()).replace(".py","_"+dateString+"_Debug.txt")
stopfile = inspect.getfile(inspect.currentframe()).replace(".py",".stop")
newdebugfile = False

codeDir = os.path.dirname(inspect.getfile(inspect.currentframe())).replace("\","/")

def printit(inText):
    global newdebugfile

    print inText

    if os.path.exists(debugfile):
        if (newdebugfile == False):
            tmpfile = open(debugfile,"w")
            newdebugfile = True
        else:
            tmpfile = open(debugfile,"a")
    else:
        tmpfile = open(debugfile,"w")

    tmpfile.write(inText)
    tmpfile.write("n")
    tmpfile.close()
    newdebugfile = True

stopFileExists = False
printit("Code Directory: "+codeDir)
printit("Starting at: "+datetime.date.today().strftime("%Y-%m-%d_")+datetime.datetime.now().strftime("%H:%M:%S"))
printit("Stopfile : "+stopfile+"/n")
while (stopFileExists == False):
    for iFile in glob.glob(codeDir+"/test_*.py"):

        thisStartfile = iFile.replace(".py",".start")

        if not (os.path.exists(thisStartfile)):
            printit ("Launching: "+iFile)
            iTmpfile = open(thisStartfile,"w")
            iTmpfile.write("started")
            iTmpfile.close()
            os.system("Start "+iFile)

    if (os.path.exists(stopfile)):
        stopFileExists = True
    else:
        time.sleep(60)

    printit("nEnd of Loop: "+datetime.date.today().strftime("%Y-%m-%d_")+datetime.datetime.now().strftime("%H:%M:%S")+"n")    

printit("Done!")

Feature classes and Tables with names starting with “nd_”.

Random luck me to discovering a bug related to feature classes whose names start with “nd_”.  It appears that you are allowed to create feature classes starting with “nd_” but ArcCatalog will not display them.  Further research shows this behavior also occurs for table and for ArcSDE (PostGres) geodatabases,  personal geodatabase, and file geodatabases–I am using ArcCatalog 10.0.

I first noticed something odd was occurring while importing a series of shapefiles into a geodatabases.  After importing 15 shapefiles, I only had thirteen feature classes despite receiving no errors during the process.  The two shapefiles that failed to import were named ND_oil_and_gas.shp and ND_Bendix_Study.shp.  Subsequent attempts to import them individually returned an error “Invalid Target Name”.

I discovered in pgAdmin III (Postgres SDE Geodatabase) that the table existed and there was an entry in sde.sde_layers for the feature class but ArcCatalog refused to show it.

I used some un-supported methods to try to resolve the problem and despite some sweating, I failed to find a way to get ArcCatalog to display these feature classes.  I did, however, at least found a way to delete them–arcpy can detect that the feature classes exists so it is able to delete them.

At least by deleting them, I can prevent leaving “invisible” feature classes from hanging out in my geodatabase.

I suspect the problems stems from how ESRI has implemented the Network dataset table-naming structure –dirty areas are stored in tables named nd_<itemid>_dirtyareas  and nd_<itemid>_dirtyobjects.  Possibly the developer  working on the ArcCatalog GUI ended up suppressing showing feature classes and tables whose names start with “nd_”.

And, just for posterity’s sake, here is a python code snippet listing the feature classes in a workspace:

import arcpy

arcpy.env.workspace = “c:/temp/_nd/F.gdb”

print arcpy.env.workspace
for fc in arcpy.ListFeatureClasses():
print fc

print “Done!”

Renaming Raster Dataset and arcpy.Exists()

Discovered something today. I was working on an arcpy script that copies a raster dataset from a file geodatabase into a Postgres SDE geodatabase and then does some boring routine tasks–building stats, creating a mosaic dataset, adding the raster to the mosaic dataset and making a couple referenced mosaic datasets.

It sometimes has trouble with the initial step of uploading the raster because of the sheer size of if (1m elevation raster for counties) and it failed today on one. It failed today so I used the ArcCatalog GUI to copy the raster and renamed it.

I then proceeded to run launch my script. Before each step, I use arcpy.Exists() extensively to check to see if various items exist before I attempt to create them. It was continuously reporting that my raster set did not exist even though I could see it in ArcCatalog.

Finally, I realized that I needed to close ArcCatalog before arcpy recognized the fact I had renamed something. To note, I was running arcpy from a separate PythonWin window, not from the ArcCatalog session I had renamed the raster dataset with.

Once I closed ArcCatalog, arcpy recognized the renaming and life was good.

I’m also suspicious now about a problem I often have running statistics on my rasters.  The ArcTool reports no errors when I create them but for some reason the raster does not show that it has statistics afterwards.  I normally have multiple ArcApplication sessions open and now suspect that perhaps this problem is due to sessions not letting go of the connection.  Stay tuned for further developments on this.

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

Change Detector arcpy Script

During a process I was working on, I needed to compare a feature class before and after some edits.  I did not quickly find anything in ArcToolbox but searching ArcResources led me to Change Detector script by Bruce Harold.  After making a couple of tweaks–for some reason in one of my feature classes, the Shape field had an upper case “S” and in the other it was a lower case “s”.  I also discovered that it needs to export to the same format (personal geodatabase, file geodatabase, shapefile) as the source data (or at least one that uses the same field name deliminator).

After minor adjustments, though, it worked like a charm.  I’ll be submitting the changes I made to Bruce and let him incorporate the changes into the official code.

FOLLOW-UP: Mr. Harold quickly responded to my email & made the change (although I haven’t checked it). Way to go Bruce!  Thanks for a handy script.