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.

Advertisements

Checking to see if a Field Index Exists Using Arcpy (ArGIS 10.0) redux

I’ve previously posted python code to check if a field index exists for both ArcGIs 9.3 and ArcGIS 10.0.

Recently I have been working on a process that was using this code but it was not working because it looks for an index with a specific name.  It was not working in this case because the name of the indexes was getting incremented as they were being created.  For example, I was building an index on the table C5ST, field RelateId ([C5IX].[Relateid]) named I_C5IX_RelateId.  That worked fine until we switched our process so now we keep multiple versions of some tables, each with a date-based suffix.

We now have tables name C5St_20110625 and C5St_20110626–the Index-name scheme, however was still creating I_C5IX_RelateId and it worked great on the first one.  But when it created the second one, even on a different table, it was automatically name I_C5IX_RelateId_2 even though the name I_C5IX_RelateId was used when trying to create the index.

Before generating relates, our code checks to see if the key fields are indexed, and if they are not, builds  an index.  Because of the naming situation, multiple, duplicate indexes were being created.  Probably not too harmful but it is a little messy.

So I re-wrote the code so that you pass the function the table name and field name that you want to check and it checks to see if there is an index existing for that field and return a Boolean.  The one little wrinkle I put in is to account for indexes that span multiple fields–the ” if (iIndex.fields[0].Name.upper() == fieldname.upper()):” statement is checking the index to see if it is on a single field or multiple fields.

 

def fieldHasIndex(tablename,fieldname):
if not arcpy.Exists(tablename):
return False

tabledescription = arcpy.Describe(tablename)

for iIndex in tabledescription.indexes:
if (len(iIndex.fields)==1):
if (iIndex.fields[0].Name.upper() == fieldname.upper()):
return True

return False

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.

Using arcpy to List Domains Assigned to Featureclass Fields

I was making an edit (adding leading “0”s) to a coded-value domain in an SDE database and realized that my edits were changing the order of the rows of my domain.  Rows were moved to the bottom of the list when they were edited.

So I went through the process of converting my domain back to a table, made my edits in Access and exported the rows to a .dbf in the order I wanted them.

I removed the domain from the field I knew it was assigned to but when I tried to delete the domain, I received an error (The domain is used as a default domain).

The Google Machine led me to an ArcForums post by  with some python code for listing all the fields with a domain.

I tweaked the original a bit, first because it was only examining feature classes in a feature dataset, not stand-alone feature datasets.  And secondly, I updated it to use arcpy.  I posted both the 9.3 version and the 10.0 version, but here is a quick look.  The key addition is the ‘listfc(“”)’ line that is the first line of the def listds() module.

import arcpy

min_workspace = "C:\\Users\\mranter\\AppData\\Roaming\\ESRI\\Desktop10.0\\ArcCatalog\\min.minstaff.sde"
#min_workspace = "C:\\TEMP\\kurst-10.mdb"

infgdb=(min_workspace)
arcpy.env.workspace = infgdb

def listfc(inDataset):
featureclasses = arcpy.ListFeatureClasses("","",inDataset)
for f in featureclasses:
print "feature class: ",f

lfields=arcpy.ListFields(f)

for lf in lfields:
if lf.domain<>"":
print "      domain",f, lf.name, lf.domain

def listds():

listfc("")

datasets=arcpy.ListDatasets ("","")
for d in datasets:
print "  dataset: ",d

listfc(d)

listds()

Zipping a File Geodatabase using Python

Ever since the ever-popular post, Zipping a shapefile using python, came out, people have been asking (one person, yesterday) for a sample of how to zip a file geodatabase using python.

The key trick, as shown in line 17, is appending the basename of the file geodatabase (“nfg.gdb/” in my example) in front of each file as you write it to the zipfile.

UPDATE: WordPress messes with the spacing when I post code, making it difficult to post code that can just be copied & pasted and have work.  So I have posted a the code HERE for downloading.

import os
import zipfile
import glob

infile = "c:/temp/nfg.gdb"
outfile = "c:/temp/nfg.zip"
def zipFileGeodatabase(inFileGeodatabase, newZipFN):
if not (os.path.exists(inFileGeodatabase)):
return False

if (os.path.exists(newZipFN)):
os.remove(newZipFN)

zipobj = zipfile.ZipFile(newZipFN,'w')

for infile in glob.glob(inFileGeodatabase+"/*"):
zipobj.write(infile, os.path.basename(inFileGeodatabase)+"/"+os.path.basename(infile), zipfile.ZIP_DEFLATED)
print ("Zipping: "+infile)

zipobj.close()

return True

zipFileGeodatabase(infile,outfile)