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

Referenced Mosaic Datasets Permissions on SDE

While I’m enjoying the functionality of the new Referenced Mosaic Dataset that have been introduced in ArcGIS10, something that I stumbled over recently was administering the privileges to a referenced mosaic dataset.

A referenced mosaic dataset is a raster datatype that uses a mosaic dataset as a base.  For example, we loaded a series of DEMs into a mosaic dataset and then created a referenced mosaic dataset that does the meters to feet conversion for us.  This saved us from having to process each DEM before loading it.  We also created a referenced mosaic dataset hillshade from the original meter data.

The problem I came across, and this was more my oversight than anything, but I had a heck of a time trying to grant privileges to the referenced mosaic dataset.  The solution was actually very simple–I needed to also grant access to the mosaic dataset that was being used as the base.

That makes some sense, you should not be able to grant access to a dataset (the referenced mosaic dataset) that would inadvertently be exposing data (the base mosaic dataset) that a user does not have access to.  This, however, conflicts a bit in that the permissions on individual rasters that make up a mosaic raster dataset are ignored.

ArcMap Field Calculator: Calculating Geometry using arcpy

One of the things I had not gotten around to doing in ArcGIS 10 was figure out how to directly manipulate the geometry of a record using the Field Calculator.  When I stumbled upon a bug in the way the Extract Values to Points tool handles Null geometries, I figured it was time to figure it out.  Setting the X, Y to 0,0 was sufficient for my needs.

I set the Parser to Python and the formula was simple once I figured out the syntax:

pPoint = arcpy.Point(0,0)

Then, to extend my knowledge, I wanted to see how to calculate the geometry based off of two fields (X and Y).   The only real trick is knowing the bracket field names using exclamation points instead of brackets:

arcpy.Point(!X!,!Y!)

So turns out everything it pretty easy and straight-forward.  Whew!

Extract Values to Points (Spatial Analyst) Bug

One of the Spatial Analyst tools we often use in ArcGIS is the “Extract Values to Points” tool.  This allows us to take a point file (well locations in our case) and attach a value (elevations) from a raster image (a DEM) to each point.

Today I was running it for the first time against an Image Service we recently published and I received a warning message,”WARNING 000957: Skipping feature(s) because of NULL or EMPTY geometry”.  But the script seemed to run and the final results said “Succeeded” so I thought it was probably fine.  But as I double-checked, I realized the results were wonky.

Turns out that I had two records with Null geometry in my point file of 397 records.  These two records threw the above error but actually had a value in the [Rastervalu] field.  Turns out all 397 records had values.  These two records were consecutive–let’s say the 100th and 101st records in my shapefile.  What happened is record 100 got the value for record 102, record 101 get the value for 103, record 102 (which has valid geometry) had the value for 104.  This pattern, each record having the value for the record 2 place after it, continued until record 396 which had the value for record 397.  Record 397 also had the value for 397.  So the final three records all had the value for the final record.

What I would have expected would be for the two records with Null Geometry to have null values in the [Rastervalu] field and the rest of the records to have the correct values.  Despite the warning, it is very misleading for all the records to end up with a value.

I have a simplified example below.  I made a point shapefile with four records.  The first, third, and fourth  records have valid geometries; the second has Null geometry.  The second record ends up with the value for the third record.  The third record, has the value for the fourth.  The fourth record being the last record, ends up with the last valid value, which was its own.

The results that I would have hoped for would be for the third record to have a Null value.

The way I envision what is occurring behind the scenes is this:  the process makes a list (more of a stack in programming terms) of result values as it processed the points but just assumes that every record will return a value so it does not track which value goes with which shape.

When it reached the two null geometries, it threw an error but continued on.  It did not add a value for these records to the stack of values–when it comes across records with valid geometry but do not intersect the raster it adds a psuedo-null value of (-9999) to the stack.  After it processed all the records it had 395 values in the stack.  It then went, one-by-one through the stack and populated the records in the output shapefile, the first record got the first value in the stack, the second record got the second value, the 100th record got the 100th value (which came from the location of the 102nd record) and so on.  At the end, the final two records received the last valid value.

This final behavior–using the last valid value–corresponds a bit to a behavior we’ve seen with ArcObjects in general.  When iterating through a table, if a field is Null for a specific row, the value from the last non-Null value for that field is often returned.

I’m in the process of submitting a bug to ESRI.  I’m not sure if this existed prior to ArcGIS 10.0 (I’m guessing it did) or if it occurs in other processes (I’m guessing it does).  I did find out that the “Extract Multi Values to Points” works as expected.  I’m guessing it is because unlike the “Extract Values to Points” which creates a new shapefile, this tool appends fields to the existing shapefile and presumably processes records one-by-one without putting the results in a virtual stack.  The “Extract Multi Values to Points” tool also does not throw any warnings.


Multiple geodatabases and ArcSDE services from one PostgreSQL server.

To better organize our ArcSDE data, we wanted to create multiple geodatabases and multiple ArcSDE services using one PostgreSQL database cluster (a cluster containing 1 machine at this point).  A side question is why can’t tables and raster be placed in Feature Datasets?  This wouldn’t be an end-all solution for what we want to do but there are some messy consequences of this limitation.

ESRI has instructions on Setting up multiple geodatabases in one PostgreSQL database cluster on Windows which was helpful but we repeatedly got an “The ArcSDE Repository was unsuccessfully completed.  Would you like to view this error?” error during the Repository Setup step, step two of four, of the ArcSDE Post Installation process. The odd thing, however, was that the wise_err.log was blank so I had no clue on why it was failing.  After much Googling and head-banging, I came across a post from Kim Doty on ESRI’s forums reporting the same problem but they were able to create a service by just continuing through the process .  Figuring I had nothing to lose, I thought I might as well try to go through the entire ArcSDE Post Installation process and see what happened.  Although I received blank error messages, I was able to successfully create the second ArcSDE instance.

I will touch on some of the highlights of the process.  In my case, I was creating a database & service with these parameters:

  • Database name: mgsgdb_lidar
  • Service name:  mgs1_sde
  • Tablespace folder: D:mgsdb1dmgsgdb_lidar

Following ESRI’s instructions, I first made custom copies of dbinit.sde, dbtune.sde, and giomgr.defs that contain the service name (mgs1_sde) within the file name.  These copies were named dbinit_mgs1_sde.sde, dbtune_mgs1_sde.sde, and giomgr_mgs1_sde.defs.

When I have trouble with the ArcSDE Post Installation, I like to go through the four steps individually so I select a “Custom” install.

The four steps making up the ArcSDE Post Installation process are shown on the next dialog.

I was able to complete the first step, Define SDE User Environment, without any problems.

Make sure to set your database name, default tablespace, and tablespace folders are distinct from your first installation.

During the second step of the ArcSDE Post Installation process, Repository Setup, make sure to use the custom files you made earlier when you come to the ArcSDE configuratino files dialog.

Later in the Repository Setup, make sure to use the correct database name:

This is the first error message I received in the process (still during Repository Setup):

Clicking “Yes”, however, shows an empty wise_err.log file:

After viewing this, I canceled out of the Repository Setup step.

Taking a look in pgAdmin III, however, I can see that both my database (mgsgdb_lidar) and tablespace (mgsgdb_lidar) were created:

I received the same error message during the third step, Authorize ArcSDE, as in the second step:

But running the fourth step, Create ArcSDE Service, resulted in this sweet message:

After editing my pg_hba.conf and opening the port in my firewall, the service was created, running, and visible.  So far, it seems to be fully functional and without problems.

ArcGIS 9.3.1 to ArcSDE 10 Connection Error

We finally installed an instance of ArcSDE 10 today.  My first attempt at connecting in ArcCatalog 9.3.1 failed with the following error:

Failed to connect to the specified server.

This release of the GeoDatabase is either invalid or out of date. [Please run the

ArcSDE setup utility using the -o install option.]

DBMS table not found[sde.sde.GDB_Release]

 

Turns out the solution was simple, this article points out that Service pack 2 is required.  After updating my ArcGIS, the problem was solved:

I was the proud, new recipient of a non-misleading error message.  After some minor head-slapping, I realized that 9.3.1 and older versions of ArcGIS are not able to connect to ArcSDE 10.  Doh!  My bad–I should have read ESRI’s clear information on this topic.  At least I never actually tried to install ArcSDe with the -o option as the first error message instructed me to do.

Displaying Node Dangles in Arcedit

Since the name of the blog is Node Dangles, I get several hits daily from searches on “Node Dangles” and I have no  information on node dangles.  This post is the first in a series to change that.

First let us, by us, I mean ESRI, define what a node dangle is.  Their online glossary actually defines a dangling arc, a dangling node is a node (endpoint of an arc) that does the non-connecting mentioned below:

“An arc having the same polygon on both its left and right sides and having at least one node that does not connect to any other arc. It often occurs where a polygon does not close properly, where arcs do not connect properly (an undershoot), or where an arc was digitized past its intersection with another arc (an overshoot). A dangling arc is not always an error; for example, it can represent a cul-de-sac in a street network.”

Matt’s definition is “a node that isn’t connected to any other node”.  The Arcedit screen shot below shows a coverage of roads and its arcs (white lines), nodes (white squares) at the ends of the arcs, and node dangles (red squares) indicating what nodes are dangling.  Too bad clearly flagging dangling chads isn’t so eary .

Displaying (and illustrating) how to display dangling nodes in Arcedit is simple, the command is “drawenvironment node dangle” or, as I’ve done above, “drawe node dangle”.  To display the dangles in a different color than the default white, follow-up with the command “nodecolor dangle 2”.  The 2 is ESRI’s value for red.

Zipping a shapefile using python

UPDATE:

After receiving a request to modify the code to ignore .lock files, I have an updated to this post.

One of the tasks I’ve been automating is publishing a weekly data update to a website. The update consists of shapefile. The trouble with shapefiles is they consist of 3 or more files with the same basename but different extensions in the same directory.

Not an overly complicated situation but a common one that ArcGIS does not have a solution out-of-the-box. Below is a bare-bones code snippet to do it. It has both the input shapefile and output zip file hard-coded. The call to the subroutine that does the work is: zipShapefile(wellsShapeFile,wellsZipFile)

import zipfile
import sys
import os
import glob
wellsShapeFile = "C:/cwi5_bk/WELLS/wells.SHP"
wellsZipFile = "C:/cwi5_bk/WELLS/test5.zip"

def zipShapefile(inShapefile, newZipFN):
   print 'Starting to Zip '+inShapefile+' to '+newZipFN

if not (os.path.exists(inShapefile)):
   print inShapefile + ' Does Not Exist'
   return False

if (os.path.exists(newZipFN)):
   print 'Deleting '+newZipFN
   os.remove(newZipFN)

if (os.path.exists(newZipFN)):
   print 'Unable to Delete'+newZipFN
   return False

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

for infile in glob.glob( inShapefile.lower().replace(".shp",".*")):
   print infile
   zipobj.write(infile,os.path.basename(infile),zipfile.ZIP_DEFLATED)

zipobj.close()
return True

zipShapefile(wellsShapeFile,wellsZipFile)
print "done!"

Another TopoToRaster Error

Subtitled:  Why error messages are good.

Came up with another error while running TopoToRaster but this time ArcGIS gave an error message that led to a solution.  Turned out all my contour lines had an elevation of 16 which TopoToRaster did not like.  I had intended to increase the elevation and inadvertently set them all to sixteen.  I had saved the previous values before editing so it turned out to be a simple fix and I didn’t have to spend a day trying figure out what was wrong.


TopoToRaster Error

In running an automated process, I had a TopoToRaster repeatedly fail on me. The only input theme was a contour theme. The process ran fine when I used the envelope of the contour theme as the output extent but when I changed it to the envelope of a polygon theme, it would bomb. The polygon’s envelope was smaller than the contour theme.

ArcCatalog would bomb out without presenting any sort of useful message.

I determined eventually that If I set the left extent as 435210, the process would work if set the right extent to 446655 or greater but would bomb if I used 446654 or less. If the the left extent was 435209, the right extent had to be 446654 or greater.

The image below shows the contours (with vertexes) and the polygon theme I was using. The two circle graphics show the approximate extent I was clipping to.

I eventually wrote python script to do a series of tests with the intent of determining the pattern of what does and does not work. The unexpected benefit was I actually got some sort of error message back. Turns out it had to do with my lack of vertices–I densified my polylines and it ran fine. Something I should have thought of earlier but it was difficult lacking any feedback.