Quick & Dirty arcpy: Compare Feature Class Table Schemas

I’m in the process of rewriting a process, moving most of the processing from arcpy to postgresql-enabled python (love me some psycopg2).

One of the QC checks I’m doing at the end of this re-write is just verifying that the feature class schemas are the same (or that the differences are intended)  under the new process as they were in the old process.

And while ArcGIS does have a good tool for this, there were a couple tweaks I wanted to make. Most notably, I wanted a list of fields that are not in both feature classes.

ArcGIS Table Compare

So I made a quick & dirty script to do that, nothing especially clever but I’ve found it useful. Download it from GitHub. I have it currently set up to work on feature layers but you should be able to change the toolbox parameter types to allow feature classes or tables.

import arcpy,sys,os

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

featureclass1 = sys.argv[1]
featureclass2 = sys.argv[2]

tableheaders = 'name, type, width, precision, domain'

def makeFieldDict(inFC):
    d = arcpy.Describe(inFC)
    printit("Dataset: "+d.baseName)
    printit("Type: "+d.dataType)
    printit("Path: "+d.catalogPath)
    printit(" ")

    lFields=arcpy.ListFields(inFC)

    printit (tableheaders)
    fieldDict = dict()
    printit (lFields)
    for lf in lFields:
        fieldDict[lf.name] = [lf.name,lf.type,lf.length,lf.precision,lf.domain]
        printit (lf.name+", "+lf.type +", "+str(lf.length)+", "+str(lf.precision)+", "+lf.domain)
    return fieldDict

fieldDict1 = makeFieldDict(featureclass1)
fieldDict2 = makeFieldDict(featureclass2)
errorList = []
printit(" ")
printit(" ")
printit("Comparing Fields:")
for iField in sorted(list(set(fieldDict1.keys()+fieldDict2.keys()))):
    if not (fieldDict1.has_key(iField)):
        theResult = " {0} not found in {1}".format(iField,featureclass1)
        errorList.append(theResult)
    elif not (fieldDict2.has_key(iField)):
        theResult = " {0} not found in {1}".format(iField,featureclass2)
        errorList.append(theResult)
    else:
        if (fieldDict1[iField] == fieldDict2[iField]):
            theResult = " {0} OK".format(iField)
        else:
            theResult = " {0} Have Different Definitions \n   {1}: {2}\n   {3}: {4}".format(iField,featureclass1,fieldDict1[iField],featureclass2,fieldDict2[iField])
            errorList.append(theResult)

    printit( theResult )

printit(" ")
printit(" ")
if len(errorList) == 0:
    printit("GOOD! No difference Found!")
else:
    printit("These Differences Found:")
    for iError in errorList:
        printit(iError)

printit("Done!")

Friday Fave: Geodatabase Geek

This Friday Fave is more for utility than pleasure.

Unfortunately, I have been working to determine why my views and query layers perform so much worse than directly accessing my feature class.

My Googling led me to Geodatabase Geek, by Trevor Hart, Eagle Technology Group Ltd.  Trevor has some real good information about Geodatabases and also  gave a good lightening talk on Usage Reporting on ArcGIS 10.1 for Server at the 2013 ESRI International Developer’s Conference.

One tool he pointed out was Mxdperfstat for benchmarking the performance of your MXD. Trevor used it to compare the performance of a Feature Class vs Query Layer vs Spatial View. While the official version is available for ArcGIS 9.3 through 10.2, I do want to point out Hussein Nasser’s 10.1 version which he put out before the official 10.1 version came out (it’s not really a version, more of a work-around but I like his ingenuity).

My results were significantly different on our 10.0 database server, the spatial view I was testing was much slower.  The query for both the spatial view and query layer was simply “Select * from featureclass

So not sure what to make of the performance yet, I’ve got a spatial index made so not sure what else I can try.

ArcSDE 10.0 Performance
ArcSDE 10.0 Performance

 

Unsupported Arc: “Rebox”ing or updating the extent of a feature class.

I’ve found that sometimes I can not find the answer to a question until I know the answer & then it becomes ridiculously easy to find the answer.

One small annoying thing that I never spent much time was when you delete features from a feature class making it significantly smaller but the envelope does not get re-sized so the zoom extent (still the original extent) is too large. This often happens to use when we convert tables to an XY theme and there are blank records–most of our data shows in Minnesota but there are some in Oklahoma (I think). Once we eliminate or correct the blank records, our data view still pops out to include a large section of the United States even though we only have data in Minnesota.

A long, long time ago, Workstation ArcInfo had a simple command, Rebox, for just this purpose (actually it still does, I just don’t get to use it anymore)–it shrunk the extent to the smallest rectangle required to enclose all the data. Up until today, I thought the request for this feature was completely ignored.

While researching something else, I was digging around in the sde tables and found one, sde.sde_layers, that had the interesting fields, minx, miny, maxx, and maxy. My quick & dangerous test (I performed it on a throw-away feature class in a throw-away geodatabase) gave me the results I wanted–once I loaded the feature class into ArcMap, the extent was a nice, tight rectangle around my features.

Is this a supported way to Rebox the extent? No.

Is it recommend by ESRI or me? No.

Will it screw up your entire geodatabase, making you lose all your data & costing you your job? Probably not but do you want to take that chance?

Will it get the job done? Maybe.  But in the process of writing this post, I found two safer ways to go about it. First, the straight-forward, sde command-line way that probably always existed that I never found until today, sdelayer -o alter had an -E option to reset the extent, including the ability to either specify it or have sde calculate it. Ok, that is usable for one person in our organization.

Previously, we had found either a VBA or other tool for doing this but had minimal success with it. Today, I found an ArcGIS 10 Add-In that is suppose to do the same thing. In my experiments (sample size n=1) it worked perfectly. If you need this sort of functionality, I would recommend trying out this Add-In first, if that fails go the sde command line route. Use the direct SQL method at your own risk!

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

Building a local, permanent cache of Bing Imagery.

I recently had an internal request to capture and store the Bing imagery for an area for future use. The user was interested in some specific images that were taken after a fire, making the ground surface-and certain geological features-much more visible. His concern was that in the future this imagery might get replaced with updated imagery taken when the vegetation has grown back.

Since it is unknown when/how this data might be used by us, we mostly wanted to capture it now & find a way to use it.

While we initiated the process of finding out what agency the data was available through, we also came up with a quick & dirty way to download the data.

Since ArcGIS 10 has made the process of loading cached basemap data a trivial process through ArcGIS Online, I have not used it much since taking a quick first look at it in 2010.

After removing my old, forgotten version and installing the latest, shiniest version of ArcBruTile, I verified it was able to display the imagery we wanted. ArcBruTile can be used to “display maps from OpenStreetMap, Bing, SpatialCloud, MapQuest, Europa Technologies, VR-TheWorld Online, Mapbox, Stamen Design and others in ArcGIS Desktop”. The cool thing for me was it builds a local cache in an open format–a bunch of jpeg files in a directory structure. All I had to do was clear the cache, and pan through the area of interest at the desired scale.

I could either spend many long boring hours manually panning, go through the process of renting a chimp to do it for me, or write some code to do it for me. I ended up making a fishnet of the area of interest and wrote a python script to pan through the area (to be posted).

After I had the images, I ended up build a Mosaic Dataset and added the images to it.  The last trick I that I had to figure out–and really I just found in it ArcForums–was how to create a mosaic dataset using relative paths. Can not be done, at least in 10, but by using the “Repair…” option to reset paths, you can make the mosaic dataset portable enough that if the reason you wanted to use relative paths was so you could move the data around or to other machines, you can. Just need to repair the paths.

So now, until we can actually track down the original data, we at least have a usable, archive of the imagery we wanted to preserve and have a way to access it in the field in a non-connected environment.

ArcMap Field Calculator: Identifying Unique Cases, Single Field

Seems like a lot of people are finding the ArcMap Field Calculator examples that I have posted useful so I will make an effort to post more of them. Most posts are generated after I do something and think that others might want to know how to do it. (Or so I can go back and remember how I did something without re-inventing it).

Something I did today was create a field (!Case!) and then populated this with a unique identifier for each different value (case) that occurred in a different field (!Feature!).

Note: python’s index statement is a 0-based search so the first case will have the value 0, the second will have 1, and so on. If you want to start the results at 1, you can make the last line: “return caseList.index(inValue) + 1”.

The basic structure for this is shown:

caseList = [ ]

def returnCase(inValue):
   global caseList

   if not inValue in caseList:
      caseList.append(inValue)

   return caseList.index(inValue)

ArcMap Field Calculator: Beware of Integer Division!

Apparently, if you post one time about ArcMap field calculator, you’re bound to get additional questions.  After my recent post about using field calculator to convert text values to numeric, someone asked about a problem they were having with another calculation they were having.

The underlying problem was that python 2.6, which is installed with ArcGIS 10, uses integer division when both the numerator and denominator are integers. The result of integer division is an integer rounded towards negative infinity.

If you’re not aware of this (or forget) then you open yourself up to some unexpected results–they can be especially hard to catch if you are using it within a larger block of code.

In this example, I’m calculating a percentage but the result is 0 for all the records because of the rounding.

The easiest thing to do in this simple example is just convert one of the two value to a non-integer value. This can be done by multiplying by 1.0 (not just “1”, you need to include the “.0”) which is a float datatype. Multiplying one of the values by a float makes that value a float and integer division no longer applies & we end up with a happy GISer.

Another option if you are using a Codeblock is to include “from __future__ import division” in your code block. Python is slowly moving away from using integer division and the ___future___ module overrides the default behavior.

ArcMap Field Calculator: Text to Double

Received a request yesterday asking how to use the ArcMap Calculator to copy values from a Text field to a Double field using python syntax.  As any good blogger would do, I immediately thought, “Awesome! Someone’s question is the perfect topic for a new blog post”.

The python parser is actually pretty good at casting values on the fly so if the values in your text field (!Day! in my example) are valid values that can be converted to a Double value, it is as simple as just setting the formula to be the text field. In my example case, I wanted to copy the value from !Day! to !DecDay! so I set the formula to be DecDay = !Day!.

That should work fine if you have clean values in your text field. In the example above, you might notice I had a selected set of 3 records that all had numeric values in the !Day! field. When I included the fourth row, which does not have a numeric value in the text field, I get this error message (“There was a failure during processing, check the Geoprocessing Results window for details.” when I use the same formula. Time to add in an error exception.

For more advanced logic, the Field Calculator dialog allows you to use a python function if you check on the “Show Codeblock” option.  In the “Pre-Logic Script Code” area (Seriously, who at ESRI came up with that name?) I entered the following function. If the value in my text field (!Day!) can be cast to a number of type float, that value is returned. If the cast is unsuccessful (IE the value in !Day! is not a number), then I return -99.

def toNum(inValue):
   try:
      outValue = float(inValue)
      return outValue
   except:
      return -99

Then in the formula portion of the dialog, I call the function, passing the value in the !Day! field: DecDay = toNum(!Day!).

Now, if you would prefer not to set all the records with non-numeric values to be -99 or other error value, not return anything. To do this, I replaced the “return -99” in the original function with a filler line (“doNothing = 4”) since the try block needs an non-empty except clause.

def toNum(inValue):
   try:
      outValue = float(inValue)
      return outValue
   except:
      doNothing = 4


And that should leave the values in the double field unscathed in your records with non-numeric values in the text field.

Shameless Plug: Check out my other blog posts on using ArcMap’s Field Calculator to calculate geometry and converting a date value to an 8 digit numeric value.

Domain Sorter Add-In Version 1.1

Almost a year ago, I updated ERSI’s Domain Sort code for VB 6 to work with ArcGIS 10. Recently, I had a comment that this Add-In caused ArcCatalog to explode if you had an open OLE connection. When I tested it, it turned out the reports were accurate.

I got around to adding in a Try-Catch around the offending chunk of code & it is now better than ever. You can download just the Add-In or the Add-In with source code or get it from ESRI’s ArcGIS Resource Center.

ArcMap Field Calculator: ArcPy Date to Decimal Function

One of the standards in our databases is to store dates as 8-digit integer values in the format of yyyymmdd. This requires us to occasionally convert values from date fields into this format.

We can do this in the ArcMap Field Calculator using this arcpy function:

def datetodouble(inNum):
     splitList = str(inNum).split("/")
     return  splitList [2]  +("0"+ splitList [0])[-2:]  +("0"+ splitList [1])[-2:]