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()

 

arcpy.ListFeatureClasses() Bug

I was recently re-evaluating our back-up procedures and discovered and found a nasty bug with the arcpy’s ListFeatureClasses request. If you have a feature class in a feature dataset with the same name, ListFeatureClasses may not find it or anything else in that feature dataset.

Unfortunately, we recently made our daily backup a python-based system that uses ListFeatureClasses and got bit by this bug.

After discovering missing data in our backups, I reconstructed what happened and found this bug. Below is arcpy code that iterates through the feature datasets in a geodatabase and lists the feature classes:

import arcpy

def copyAll():
    for iFeatureClass in arcpy.ListFeatureClasses():
        print(" Feature Class: {0}".format(iFeatureClass))
    iFeatureClassFull = None

testGDBname = "mgs_sandbox.sde"
arcpy.env.workspace = testGDBname

copyAll()
for iFD in arcpy.ListDatasets("","Feature"):
    print("Feature Dataset {0}:".format(iFD))
    arcpy.env.workspace = testGDBname+"/"+str(iFD)
    copyAll()

And here is a screen shot of the contents of a test enterprise geodatabase, you’ll see it has a feature data set named “outcrops” that has a feature class also named “outcrops” within it:
sandbox

And the results list only the feature dataset:

results

But if I rename the feature dataset (e.g. outcrop_fd), the results are what I would hope for:

results results2

I found that the feature class does not even need to be within the feature dataset and also the problem does not always occur, I have had the code successfully run in some cases.

Once I confirmed the problem, I did find this thread from almost three years ago that mentions the bug. One poster indicated the same thing occurs in ArcObjects which leads me to think something may not be getting registered right in the sde tables.

I was not able to re-create this using either personal or file geodatabases.

So I adopting the policy of not using the same name for a feature dataset as for a feature class.

 

Arcpy: Check if a field exists

I was helping a co-worker who needed to check if a field exists in their arcpy script. Since we were located at their computer, I thought I would just do a quick Google search and pull the code off this blog. Seemed logical since I the original purpose was exactly that—to serve as a handy, public place to store code snippets that I use & that others might find handy.

Anyhow, my Google Search on “Node Dangles field Exists” came up with a 9.3 script to check if field index exists. And I also have a 10.0 version but did not come up with the field exists snippet. So here it is:

image

def fieldExists(inFeatureClass, inFieldName):
   fieldList = arcpy.ListFields(inFeatureClass)
   for iField in fieldList:
      if iField.name.lower() == inFieldName.lower():
         return True
   return False

ArcMap Field Calculator: Number of parts in multi-part feature

In the last week, I have looked for multi-part features a couple of times. Today, I was looking for multi-part polygons after dealing with the fall-out of a case of Clip Gone Wild as shown below.

image

I have not found a way to write a query to find these but Field Calculator does allow you to calculate a field’s value to the number of parts.

Using the Python parser, just write the formula (note that case matters): !shape!.partCount

image

Quick & Dirty python: Converting a text file to audio (.wav)

This is a bit of a tangent but for some crazy reason, I wanted to convert some text to audio so I could listen to it while I drive. A quick Google search left me without any freeware that could handle the 53 page document–there are some cool websites that do text to mp3 like vozme and YAKiToMe! but they didn’t convert the whole document. I then found pyTTS, a python package that serves as a wrapper to the Microsoft Speech API (SAPI) , which has been in version 5 since 2000. But I didn’t easily find a version of pyTTS for python 2.6. So I decided to see if I could roll my own.

As it turns out, getting python to talk using SAPI is relatively easy. Reading a plain text file can be done in a few lines.

from comtypes.client import CreateObject

infile = "c:/temp/text.txt"

engine = CreateObject("SAPI.SpVoice")

f = open(infile, 'r')
theText = f.read()
f.close()

engine.speak(theText)

And it wasn’t that much more to have it write out a .wav file:

from comtypes.client import CreateObject

engine = CreateObject("SAPI.SpVoice")
stream = CreateObject("SAPI.SpFileStream")

infile = "c:/temp/text.txt"
outfile = "c:/temp/text4.wav"
stream.Open(outfile, SpeechLib.SSFMCreateForWrite)
engine.AudioOutputStream = stream

f = open(infile, 'r')
theText = f.read()
f.close()

engine.speak(theText)

stream.Close()

And with that chunk of code, I was able to convert my 54 page document into a 4 hour long .wav file (over 600 MB) that I used another software package to convert to .mp3 (200 MB). The voice is a bit robotic but not too bad, I just hope the content that I converted (a database specification standard) doesn’t put me to sleep while I drive.

Quick & Dirty arcpy: Autopan ArcMap using arcpy

Question: How do I get ArcMap to automatically pan through an area.

As I mentioned in a previous post, I recently had the need to have ArcMap automatically pan through a project area. My first attempt was to print a series of data-driven pages (using a fishnet polygon layer as the index) this but that did not accomplish what I needed so I switched to arcpy, which made the task simple enough. Nothing special or tricky about this code, but just did not find it anywhere else.

The one thing to note is that I have a 1 second pause between pans–this was to allow image tiles to download. You will need to adjust the delay to meet your needs. The toolbox and code can also be downloaded.

import sys,arcpy,datetime
inLayer = sys.argv[1]

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

mxd = arcpy.mapping.MapDocument("CURRENT")

arcpy.MakeFeatureLayer_management(inLayer, "indexLayer")
cur=arcpy.SearchCursor("indexLayer")

df = arcpy.mapping.ListDataFrames(mxd)[0]
newExtent = df.extent

iCount = 0
iTotal = (arcpy.GetCount_management("indexLayer").getOutput(0))

for row in cur:
    thisPoly = row.getValue("Shape")
    newExtent.XMin, newExtent.YMin = thisPoly.extent.XMin, thisPoly.extent.YMin
    newExtent.XMax, newExtent.YMax = thisPoly.extent.XMax, thisPoly.extent.YMax
    df.extent = newExtent
    time.sleep(1)
    iCount+=1
    printit("Panned to feature {0} of {1}".format(iCount,iTotal))

del row
del cur