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!

Advertisements

Quick & Dirty python: Converting a text file to audio (.wav)

This is a bit of a tangent but for some crazy reason, I wanted to convert some text to audio so I could listen to it while I drive. A quick Google search left me without any freeware that could handle the 53 page document–there are some cool websites that do text to mp3 like vozme and YAKiToMe! but they didn’t convert the whole document. I then found pyTTS, a python package that serves as a wrapper to the Microsoft Speech API (SAPI) , which has been in version 5 since 2000. But I didn’t easily find a version of pyTTS for python 2.6. So I decided to see if I could roll my own.

As it turns out, getting python to talk using SAPI is relatively easy. Reading a plain text file can be done in a few lines.

from comtypes.client import CreateObject

infile = "c:/temp/text.txt"

engine = CreateObject("SAPI.SpVoice")

f = open(infile, 'r')
theText = f.read()
f.close()

engine.speak(theText)

And it wasn’t that much more to have it write out a .wav file:

from comtypes.client import CreateObject

engine = CreateObject("SAPI.SpVoice")
stream = CreateObject("SAPI.SpFileStream")

infile = "c:/temp/text.txt"
outfile = "c:/temp/text4.wav"
stream.Open(outfile, SpeechLib.SSFMCreateForWrite)
engine.AudioOutputStream = stream

f = open(infile, 'r')
theText = f.read()
f.close()

engine.speak(theText)

stream.Close()

And with that chunk of code, I was able to convert my 54 page document into a 4 hour long .wav file (over 600 MB) that I used another software package to convert to .mp3 (200 MB). The voice is a bit robotic but not too bad, I just hope the content that I converted (a database specification standard) doesn’t put me to sleep while I drive.

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: Identifying Unique Cases, Multiple Fields

You may have noticed that this post–ArcMap Field Calculator: Identifying Unique Cases, Single Field–specifies “Single Field”. Yes, that was my version of a cliff-hanger post.

The basic structure I listed in that post can be expanded on to satisfy your needs. The example in my earlier post was case sensitive for example, you could modify it so it treats “a” the same as “A”.

Today’s example groups records into different cases based off the values of two fields, !county_c! and !feature! and required only minor modifications.

The calling line was modified from:

returnCase(!feature!)

to:

returnCase(!county_c!,!feature!)

to accommodate passing both values.

The function definition likewise was modified to accept two values, this:

def returnCase(inValue1):

to:

def returnCase(inValue1, inValue2)

And this line was added, creating a list from the two values passed in:

inValue = [inValue1, inValue2]

 

(Note: The same results could have be achieved by using the original function by creating the list in the calling statement:  returnCase([!county_c!,!feature!] )

 

caseList = [ ]

def returnCase(inValue1, inValue2):
   inValue = [inValue1, inValue2]
   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.