Quick & Dirty arcpy: Bulk Changing Field Values

In mapping cross sections, our geologists often find themselves renaming their stratigraphic units midway, or at the end, of creating multiple cross sections.  This can cause a situation where we need to change multiple values in multiple fields in multiple feature classes–a situation that can get messy very fast.

Perfect situation for a quick & dirty arcpy script and, in this case, an ArcToolbox tool that can be downloaded.

This tool will change all feature classes in the O:\clay_cga\sand-distribution_model\dnrPackages\stratlines directory.

It will look at two fields, [strat] and [unit] and make these changes:

  • “go” becomes “gro”
  • “goc” becomes “grc”
  • “sgb” becomes “grb”

And since I have Case Sensitive checked, “Go” will not get changed to “gro”.  Also note that only full values that match values in the Old Value List get changed, part matches are left as-is so “got” would be left as-is even though the first two characters match “go”.

 

Bulk Field Change

import arcpy
import sys, string, arcgisscripting
import arcpy

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

def printerr(inString):
    print inString
    arcpy.AddError(inString)

def fieldExists(tablename,indexname):

 if not arcpy.Exists(tablename):
  return False

 tabledescription = arcpy.Describe(tablename)

 for iField in tabledescription.fields:
     if (iField.Name.lower() == indexname.lower()):
         return True

 return False


if len(sys.argv) > 1:
    inDirectory = sys.argv[1]
    inFieldNameRaw = sys.argv[2]
    oldValue = sys.argv[3].replace(","," ")
    newValue = sys.argv[4].replace(","," ")
    caseSensitiveRaw = sys.argv[5]
else:
    inDirectory = r"C:\temp\test\stratest"
    inFieldNameRaw = "strat"
    oldValue = "go, goc, sgb".replace(","," ")
    newValue = "gro grc grb".replace(","," ")
    caseSensitiveRaw = "true"

caseSensitive = (caseSensitiveRaw.lower() == "true")
fieldNameList = inFieldNameRaw.replace(","," ").split()

printit("Starting")
printit(" Workspace: "+str(inDirectory))
printit( " inFieldName: "+str(inFieldNameRaw))
printit( " oldValue: "+str(oldValue))
printit( " newValue: "+str(newValue))
printit( " caseSensitive: "+str(caseSensitive))

valueDict = dict()

def initialQC():
    global valueDict

    if not (arcpy.Exists(inDirectory)):
        printerr("Workspace {0} does not exist".format(inDirectory))
        return False

    if (len(oldValue.split()) <> len(newValue.split())):
        printerr("Number of values in {0} does not equal number of values in {1}".format(oldValue,newValue))
        return False

    iValueIndex = 0
    for iOldValue in oldValue.split():
        if (caseSensitive):
            thisKey = iOldValue
        else:
            thisKey = iOldValue.lower()

        if (valueDict.has_key(thisKey)):
            printerr("ERROR: Value, {0}, is repeated, cancelling...".format(thisKey))
            return False

        valueDict[thisKey] = newValue.split()[iValueIndex]
        iValueIndex+=1
    return True

def makeFieldList(inFC):
    thisFieldList = []

    for iField in fieldNameList:
        if (fieldExists(inFC,iField)):
            thisFieldList.append(iField)

    return thisFieldList


def main():
    arcpy.env.workspace = inDirectory
    printit(valueDict)
    for iFC in arcpy.ListFeatureClasses():
        printit("Working on {0}".format(iFC))

        iFieldList = makeFieldList(iFC)
        if (len(iFieldList) == 0):
            printit(" No fields to change, Skipping...")
            continue

        rows = arcpy.UpdateCursor(iFC)

        changes = 0
        printit(" Changing Rows")
        for row in rows:
            iChange = 0
            for iField in iFieldList:
                iValue = str(row.getValue(iField))
                newValue = iValue

                if valueDict.has_key(iValue):
                    newValue = valueDict[iValue]
                else:
                    if not (caseSensitive):
                        if valueDict.has_key(iValue.lower()):
                            newValue = valueDict[iValue.lower()]

                if (newValue <> iValue):
                    printit("CHANGE {0}".format(newValue))
                    row.setValue(iField,newValue)
                    iChange+=1

            if (iChange > 0):
                changes+=1
                rows.updateRow(row)
        printit(" Made {0} changes".format(changes))
        del row
        del rows

    printit("Main")

if (initialQC()==True):
    main()

printit("Done")

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.

 

ArcMap Maximum Legend Items from a ArcGIS Server service

Recently we took a call from a user who could not see the legend for one of the feature classes in one of our services. (Precambrian Bedrock in http://mgsweb2.mngs.umn.edu/arcgis/services/state/mnbdrkgeology )

After trying some standard things–restarting the service, checking the source .MXD–I turned to The Google Machine and quickly found help from ESRI: http://support.esri.com/zh-cn/knowledgebase/techarticles/detail/33741 .

Turns out the default number of legend items ArcMap will display from an ArcGIS Server map service layer is 100 and we had 102 in the problematic layer. (I plead innocence, blame the geologist for needing that many categories).

The solution was to edit the Windows registry and change the setting for “Maximum Legend Count” from 100 to something higher than 102. After doing that (see path note below) and restarting ArcMap, the legend showed for us.

Path Details
The path turned out to be a little different on Windows 7 than the paths indicated in the help article.
ESRI indicated that the path varied by ArcGIS version:

  • ArcGIS 9.3.x: HKEY_CURRENT_USER > Software > ESRI > MapServerLayer
  • ArcGIS 10.0: HKEY_CURRENT_USER > Software > ESRI > Desktop10.0 > ArcMap > Server > MapServerLayer
  • ArcGIS 10.1: HKEY_CURRENT_USER > Software > ESRI > Desktop10.1 > ArcMap > Server > MapServerLayer

I actually changed this in a number of locations–presumably once for each user profile on my machine. Each followed a pattern something like:

HKEY_USERSS-1-5-23412341234-123412342134-12341234123-1302SoftwareESRIArcMapServerMapServerLayer

Extra Information

Out of cursiousity, I wondered if the 100 item limit was a per service or per layer limitation so I set my limit at 103. Because the service has several other layers, the total number of layers in the legend was about 130. Everything drew OK so it appears the limit applies per layer.

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

ArcSDE Tip: Finding an ArcSDE instance’s port

One of my main tasks right now is to document many of the details of maintaining ArcSDE geodatabases so I anticipate having several blog posts on this topic that are re-writes of documents I am working on. I am presuming that the person will have no ArcSDE experience so I am documenting very detailed information.

Almost all of the ArcSDE commands require that you specify which instance (service/port) the command applies to by using the “-i” parameter.

ArcSDE Instance Parameter
ArcSDE Instance Parameter

Since we have multiple ArcSDE geodatabases, I like to have a handy-dandy sticky note with all the geodatabases and their respective ports on the side of my monitor.

But when that is not handy and I can never remember the ArcSDE command line syntax to get a list of instances and their ports–I mean remembering “sdeservice -o list” is difficult at my age.

sdeservicelist

The quickest and most reliable way I’ve found to get the instance number is just to check the properties of your SDE connection file in ArcCatalog, right-click on it and select “Connection Properties”.

connectionproperties

And the port is right there in the service entry (5164).

connectionproperties2

SDEINTERCEPT & SDEINTERCEPTLOC

Awhile ago, I had a ArcSDE problem that required ESRI technical support to help trouble-shoot. The problem was odd but was resolved by rebooting the server.
During the process, though, the support person had me set a couple of environment variables for logging SDE activity on the client machine.

The settings were SDEINTERCEPT and SDEINTERCEPTLOC.

From ESRI’s Help, SDEINTERCEPT specifies what activity to log and SDEINTERCEPTLOC specifies where to save the log files.

I recently deleted the directory I made for the log files but did not remove the variables and I noticed that one of my python scripts reported a weird error (but continued to run, I think). I tracked it back to these variables and realized what I had done.

Googling SDEINTERCEPTLOC did lead me to some helpful information like:

The SDEINTERCEPT blog where Ken posts ArcSDE help.

This ESRI post about troubleshotting geoprocessing problems.

And this ESRI technical article about diagnosing ArcSDE Connections.

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.