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

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

 

Friday Fave: Tapiriik

This Friday Fave is a little bit different.

My interest in geospatial technologies (although we just called it GIS back then) largely because I wanted to measure my running routes more accurately and efficiently than the paper map & scrap of paper method I was using in the early 90s. When I was introduced to GIS, I knew what I was going to use it for.

Now that GPS technology is ubiquitous–I’m currently using four different GPS devices, at the same time, on my bike rides–I seldom have to use a map to measure my routes. I may still use MapMyRun.com to plan a route ahead of time if I’m running in a new area or trying to plan a loop of a certain distance but GPS has really made it so simple to just go out and run.

While there are several GPS options available, I have used Garmin ever since their 405 Forerunner came out. This was their first watch that didn’t get confused for a Timex-Sinclair 1000 strapped to your wrist.

Timex Sinclair 1000
Timex Sinclair 1000

Garmin’s watches upload their data to Gamin Connect, which works fine, but for a project I recently started, I wanted to down load all my data which Garmin Connect does not make easy. I had over 1,000 data logs and downloading them individually was not going to happen.

A little Googling led me to tapiriik.com, which allows you to share data amongst several different online services that endurance athletes might use including Runkeeper, Strava, Garmin Connect, SportTracks, DropBox.com, and Training Peaks.

https://tapiriik.com/

 

You can use Tapiriik for free (you just need to visit their website to start synchronization) or for $2 per year they will automatically synchronize data between your accounts. You just provide your account information for whichever sites you want to synchronize and either visit their site or pay $2 and it will automatically share data between your accounts.

I linked to my Garmin Connect and Dropbox accounts and tapiriik and after some chugging, I had .GPX files for all of my data.

It was easy and didn’t take very long. Definitely met my needs–I haven’t shared data between other services but I do use the desktop version of SportTracks and considered using their web version of the software but didn’t know how I would upload all my data.  Now I know.

This product definitely saved me a bunch of time. The one thing I wonder, though, is what does “Tapiriik” means?

 

Friday Fave: Custom Maps App for Android

I admit, I love picking up freebie maps. Whether it is from the front desk of a hotel or from the bicycle shop, there is a certain appeal to seeing what people put on maps. I have maps organic orchards, breweries, Minnesota authors, rails to trails, zoos, fictional places, race maps, and a variety of other things that someone felt the need to cartographize.

http://thefriends.org/wp-content/uploads/2012/12/mn_writers_on_the_map_web_download.pdf
http://thefriends.org/

So, with all these paper maps lying around, I was thrilled to find Custom Maps, a free app on Google Play.

This app allows you to take a picture (or use an image on your device) and georeference it.  You can then view your location on the map. I’ve tried it with a few maps and have been happy with the results.

Similar to georeferencing in a desktop application, you select points on the map and then the same points on a control. You need at lest two or three points. After doing it a couple of times, the process becomes pretty easy and it just takes a minute or two from taking the picture to viewing your current location on the map.

http://www.custommapsapp.com
http://www.custommapsapp.com

An example use case scenario for this app is you are lugging your family around at an overwhelming large amusement park and never quite sure where you are. You can take a picture of the map you picked up at the entrance, georeference it, and your phone will then show exactly where you are on that map and where the nearest loo is.

The biggest limitation I’ve seen so far is that there is a limit of 3-5 megapixel size limit on the image. Apparently this is an android limitation on how much memory an app can consume. But if you adjust your camera settings not to exceed this, you should be good.

So far, I’m enjoying this app. The author, Marko Teittinen (a good Finn name), has made the source code open source so I look forward to digging into it in more detail.

 

Friday Fave: Cartozia Tales

Obviously Cartographers belong in the same category as other superheroes like Superman, Batman, and Spiderman and we finally have a our own comic book to prove it.

Cartozia Tales is a collaborative effort of nine indy artists with two guest artist each issue.

Cartozia

They have an interesting plan, they’ve split the world into nine regions (what’s the name for the ninth of an area, nona-rant?) and the artist will tell a story from a different region each issue. They may build off the region’s story from the previous issue, continue on the story they told the previous issue, and start something fresh.

After being funded by a KickStarter campaign, they’ve shipped three issues already and have at least seven more promised. Think the project definitely deserve’s a look.

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.

Sorting a Coded-Value Domain Add-In (ArcGIS 10)

I am working on an data-entry application to edit feature classes that contain several coded-value-domains. The problem with some of the domains, however, is that some entries have been added after the initial creation.  So the first 25 entries are in alphabetical order and there are some stragglers at the end that are in the order they were appended.

This can be confusing for users–they go to select “Milli Vanilli” and look between “Madonna” and “Motley Crue” but can not find their favorite band there–they have to go to the end of the list to find their selection.

In the past, I have gone through the tedious process of exporting the domain to a table, sorting the table, removing the domain from the necessary field(s), deleting the domain, re-importing the table back in a new domain and finally re-applying the domain to the necessary field(s). Let’s just say I didn’t do this until someone asked a few times and I didn’t have anything more exciting–like a root canal–I could busy myself with.

But this new application contains more domains than any of other datasets so it was time to find a better solution. ESRI does have a Domain Sort Developer Sample.  It, however, did not play nice with ArcGIS 10.

So I went ahead and update it from VB 6 to VB.Net/ArcObjects 10.  I made an Add-In that can be installed by downloading the .esriaddin file and double-clicking on it.  The source code is also available.

This will add an ArcCatalog Toolbar that can be added by going to Customize-Toolbars-Domain Sorter Toolbar.

This will add a toolbar with one button.  The button enables whenever you select a geodatabase with at least one coded-value domain.

This brings up a Windows form that lets you sort any domain by either the code or description, ascending or descending.  Once you hit “OK” it re-sorts your domain.

The only problem I have had is that only the owner of a domain is allowed to edit it on an SDE geodatabase.

But other than that, the button allows you to easily keep your domains sorted.

http://edndoc.esri.com/arcobjects/9.2/CPP_VB6_VBA_VCPP_Doc/COM_Samples_Docs/Geodatabase/Schema_Creation_and_Management/Sort_a_domain/e826c5a8-9740-4f0b-86b6-d3b834735574.htm

ArcMap Field Calculator: Calculating Geometry using arcpy

One of the things I had not gotten around to doing in ArcGIS 10 was figure out how to directly manipulate the geometry of a record using the Field Calculator.  When I stumbled upon a bug in the way the Extract Values to Points tool handles Null geometries, I figured it was time to figure it out.  Setting the X, Y to 0,0 was sufficient for my needs.

I set the Parser to Python and the formula was simple once I figured out the syntax:

pPoint = arcpy.Point(0,0)

Then, to extend my knowledge, I wanted to see how to calculate the geometry based off of two fields (X and Y).   The only real trick is knowing the bracket field names using exclamation points instead of brackets:

arcpy.Point(!X!,!Y!)

So turns out everything it pretty easy and straight-forward.  Whew!

Another TopoToRaster Error

Subtitled:  Why error messages are good.

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