Sunday, February 28, 2010

arcgisscripting breaks python print function

I've bitten the bullet and finally decided to learn how to write python scripts to work with ArcGIS. The first problem I came across was that standard output appears to be silently redirected without any warning. So when you want to do some basic poking around like:
import arcgisscripting

gp = arcgisscripting.create(9.3)

tools = gp.ListTools("*")

for tool in tools:
....print(gp.Usage(tool))
this turns out not to work, because behind the scenes somewhere the GIS processing object is playing silly-buggers with stdout, so print doesn't work like it should. You can force stdout to be 'normal' by doing something like this:
import arcgisscripting

oldstdout = sys.stdout
gp = arcgisscripting.create(9.3)

#this call to ListTools has the undocumented strange
#side-effect of changing sys.stdout
tools = gp.ListTools("*")
gp_stdout=sys.stdout

for tool in tools:
....usagestr=gp.Usage(tool)
....sys.stdout=oldstdout
....print(usagestr)

I hope this isnt a sign of what I'm in for while learning this stuff -- stange undocumented behind-the-scenes funny-business.

Sunday, February 21, 2010

importing raster data into ArcGIS

This is not so much a blog post as a reminder to myself on how to do something.

When I work with spatial (2D or 3D) data, I generally like to do my own data processing (i.e. write my own software) rather than try and write VB code to work with ArgGIS. If I just want to do something vanilla that ArcGIS handles via its point-and-click interface, I'll use that, but more than that and I prefer to process the data outside of ArcGIS.

However, it is often the case that, after doing some data processing, I'll want to import the data back into ArcGIS to view it and/or do minor tweaking of the sort that ArcGIS does well (contour lines, etc). Sometimes just getting the data back in turns out to be a pain in the arse, so I've decided to put this post up to remind myself how to do it for importing grid (raster) data.

OK, so lets suppose you have some grid/raster data (temperature, rainfall, whatever). Lets say it is in a plain text file in this format:

POINT_X,POINT,Y,ELEVATION
LONG1,LAT1,VALUE1
LONG2,LAT2,VALUE2
LONG3,LAT3,VALUE3
...


(I'm assuming a geographic coordinate system but this should work also for a projected coordinate system. Also, you can have an excel file if you want too.)

Now, we want to use this data to create a raster dataset. Here is how to do it

  1. In ArcMap, click "Add XY Data" and select the text (or excel) file with the x,y,z triples.
  2. Make sure you select the spatial reference system for the dataset (it will be 'undefined' by default)
  3. This should give you a 'point' dataset. Which may not be what you are after. Its not usually what I am after, for example. I usually want a proper raster data layer. So, we need 1 more step....
  4. In the ArcToolBox, choose "Conversion Tools->Point to Raster".
Done. This should now give you raster data that has a specified coordinate system. Other ways of importing text data to raster (ASCIIToRaster, for example) don't seem to let you specify the coordinate system, so I've ended up going through this two step process. Anyone knowing a simpler way, please let me know.