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: 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:]

ArcMap Field Calculator: Create a Unique ID

One of the common functions I have to do is assign each record in a feature class with a unique identifier–normally just a sequential number from 1 to N.  In ArcView 3.x, the formula was simply “rec + 1” if I wanted to start with the number 1.

In ArcGIS, the process got a little more complex–you had to write a little VBA in Field Calculator as described by ESRI.

While this option still exists in ArcGIS 10, I believe it will disappear when 10.1 comes out and VBA support is completely eliminated.  But it is doable using Python which will continue to be supported.

Googling around, I did not find an exact answer but Dave Verbyla, Professor of GIS/Remote Sensing at the University of Alaska has a posted some samples that served as a good starting point.

In the Pre-Logic Script Code box, I declare a variable (counter) and a function. Then in the formula, I call the function.

counter = 0
def uniqueID():
  global counter
  counter += 1
  return counter

While composing this post, I actually wanted a concatenated value; “OC” plus an 8 character numeric sequential number starting at OC00000001 so the actual code is shown below:

ArcGIS Add-In Custom Mouse Cursor

I was working on a project and wanted my own custom mouse cursor and did not easily find a way to make your own in ESRI’s instructions.  But, once you know how to do it, it is pretty easy.  In Visual Studio, Add a New Item:

Add a Cursor File:

You can edit your cursor with the editor program in Visual Studio.  Once you satisfied with how it looks, make sure that the Build Action on the cursor is “Embedded Resource”.

Then you can set your cursor with two lines of code. Not that my cursor is in my QDI.QdiAddIn Namespace:

       
Dim pCursorStream As System.IO.Stream = Me.GetType.Assembly.GetManifestResourceStream("QDI.QdiAddIn.NewCursor.cur")
MyBase.Cursor = New System.Windows.Forms.Cursor(pCursorStream)

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

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.

Quick & Dirty arcpy: Field Listings

I have to often get a table structure for a feature class or table into either a spreadsheet or word processing document.  There might be an easy way to do this in ArcGIS 10 but I haven’t found it.  So, as is my nature, I decided to roll my own.

This is a bare-bones script that iterates through the fields, printing the field name, type, width, and precision.  There are three optional features to it:

  • You can choose to have it list the domain, if there is one, on each field.
  • You can have it write to a text file (otherwise you can just copy & paste the results from the results window).
  • You can have it count the number of populated records.  This can take a long time if working with a large dataset.  Also note that my logic for determining what constitutes being populated may not be what you need but the structure is there.  I also do not account for all field types, if the field is of a type I have not account for, the code will return -999.

To use the script from ArcToolbox, you need to pass it four parameters, their Names, type, whether they are input or output, and whether they are required or optional are:

  • featureclass, Table, Input, Required
  • includedomainstring, Boolean, Input, Required (controls whether or not domains are exported)
  • doCountsRespone, Boolean, Input, Required (controls whether or not you want to get the number of populated records.  (Your definition of populated may vary from my code)
  • outputFile, File, Output, Optional (optional output file to write)

Here is the code, but you are better off just downloading it since I haven’t figured out a good way to have WordPress play nicely with python’s indenting.

# Name: ListFields-arcpy.py
#
# Purpose: Lists the fields, their type, width, and precision
# Can either have it export it to a CSV file or copy
# and paste from the results window.
#
# To use, create a tool from the script and add 3 parameters:
#  1) Table, Input, Required
#  2) Boolean, Input, Required (controls whether or not domains are exported)
#  3) Boolean, Input, Rekquired (controls whether or not you want to get the number of
#  Populated records.&nbsp; (Your defintion of populated may vary from my code)
#  4) File, Output, Optional (optional output file to write)
#
#

import arcpy,sys,os

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

if len(sys.argv) > 4:
 featureclass = sys.argv[1]
 includedomainstring = sys.argv[2]
 doCountsRespone = sys.argv[3]
 outputFile = sys.argv[4]
else:
 featureclass = "C:/temp/before.shp"
 includedomainstring = "false"
 doCountsRespone = "true"
 outputFile = "C:/temp/before.csv"

if (outputFile == ""):
 doOutputFile = False
else:
 doOutputFile = True

if (str(doCountsRespone).lower() == "true"):
 doCounts = True
else:
 doCounts = False

if (str(includedomainstring).lower() == "true"):
 includedomain = True
else:
 includedomain = False

lfields=arcpy.ListFields(featureclass)

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

tableheaders = 'name,type,width,precision'

if (doCounts == True):
 tableheaders+=",count"

if (includedomain == True):
 tableheaders+=",domain"

if (doOutputFile):
 tmpfile = open(outputFile,"w")
 tmpfile.write(tableheaders)
 tmpfile.write("n")

printit (tableheaders)
for lf in lfields:

 pThisline = lf.name+","+lf.type +","+str(lf.length)+","+str(lf.precision)

 if (doCounts == True):

 rowCount = 0

 #Note that I do not account for all field types
 #Also note that my definition of being populated may vary from yours.
 #I am using -999 as a flag to indicate a field type was not successfully
 #identified.
 if (lf.type == "Double") or (lf.type == "Single")&nbsp; or (lf.type == "Integer") or (lf.type == "SmallInteger"):
  queryString = '"'+lf.name + '" > 0'
  rows = arcpy.SearchCursor(featureclass, queryString, "", "", "")
 elif (lf.type == "String"):
  queryString = '"'+lf.name + '" <> ' + "''"
  rows = arcpy.SearchCursor(featureclass, queryString, "", "", "")
 else:
  rowCount = -999
  #rows = arcpy.SearchCursor(featureclass, "", "", "", "")

 if (rowCount == 0):
  for row in rows:
   rowCount+=1

 pThisline=pThisline+","+str(rowCount)

 if (includedomain == True):
  pThisline=pThisline+","+lf.domain

 printit (pThisline)

 if (doOutputFile):
  tmpfile.write(pThisline)
  tmpfile.write("n")

if (doOutputFile):
 tmpfile.close

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 = r"C:\Users\mranter\AppData\Roaming\ESRI\Desktop10.0\ArcCatalog\min.minstaff.sde"

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