More arcgisscripting hacking

As I mentioned in an earlier post, I've been fumbling around trying to learn how to get a python program that can talk to ARCGis. After getting around some initial irritations, I'm actually making some headway. I found this series of posts particularly helpful. As a reminder to myself, and as help for anyone else grappling with this, I'll post the first non-trivial script I've come up with here. This script is still pretty simple -- it opens up a shapefile and reads all the information associated with the shapefile. Here it is:



import arcgisscripting
import sys

#This is my first attempt to do some python
#scripting with ArcGIS Geocomputation object

#first we create the gis scripting object
gp = arcgisscripting.create(9.3)

gp.Workspace = "C:/tmp"
#allow old files to be overwritten
gp.OverWriteOutput = True

#lets print out some information about a layer
shapefile= "C:/tmp/test.shp"
desc = gp.Describe(shapefile)

#From Describe object
print desc.DataType
print desc.CatalogPath

#From FeatureClass object. The Describe object
#that is returned by the above call to Describe
#is also a 'FeatureClass' object. I havent quite
#got my head around the data model used by the
#geoprocessing object, so I don't quite understand
#how we know which fields are available and which
#arent after a call to describe.... but these work
#for the shape file I'm using.
print desc.ShapeType
print desc.FeatureType
print desc.ShapeFieldName

#list the fields. This is the same as
#returned by ListFields()
for f in desc.Fields:
....print " "+f.Name


#now go and read each line in the spatial
#data table associated with the shapefile.

#first we get a 'read' cursor
readCursor = gp.SearchCursor(shapefile)

#make sure we are at the start
readCursor.Reset()

#now go through and read (and print out) each row
row = readCursor.Next()
while (row != None):
....#print the contents of the row
....resstr = "ROW: "
....for f in desc.Fields:
........resstr = resstr + str(row.GetValue(f.Name)) + " , "
....print resstr

row = readCursor.Next()


#done!


pure eye-candy

As part of the UrbanIT research project I've been co-opted into at UNSW, I've been writing some code to 'visualize' the contents of various databases holding information about the built environment. To be very brief, the aim of the project is to integrate the existing disparate data-sets that are out there to create a single unified view of the urban environment. So, you may have data on (say) the energy efficiency of different buildings, the cost of maintaining those buildings, the price of those buildings, etc etc. All this data is tucked away in different databases and you need (essentially) an integration layer sitting on top to, well, integrate it all. Once you do this, you can do relatively complex spatial queries across multiple data sources. Convert the output to KML and you get some pretty pictures like this:



This is a mashed-up visualization of two seperate queries -- one asking for the buildings within 60m of a particular lat/long, coloured by their main use (commercial:blue, industrial:red,residential:green,other:grey). The second is a query on a particular residential building with the individual apartments coloured by their estimated value (warmer colours == more expensive apartments).

Sorry, I cant provide the KML itself, as the data associated with the buildings shown in the picture is pretty rich, and is confidential.

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.

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.

So where do people want to live?

Shown below is a picture from some work in progress, where I am looking at the location choice preferences of different household types in Sydney. The picture shows where households with young children (i.e. all less than 15 years) would like to live. The graph shows the spatial variation in preference for detached dwellings only. I do look at the effect of dwelling type, but lets do one thing at a time shall we....

You can pretty safely ignore the scale -- just think warmer colours == 'more desirable'. (If you really want to know, its the log-probability of a couple-with-kids household choosing a dwelling in that region in a discrete choice experiment.... but lets not go there).



And what does it show? Well, couples with kids like the northen harbourside suburbs (i.e. Mosman, and around middle-harbour), and big chunks of the Eastern suburbs (Vaucluse, Centennial Park, etc). No real surprises there.

Other results show that couples with young kids do not like attached dwellings -- they really do want the suburban detached home. Well, most of them anyway.. there is always a distribution of preferences.

And what about couples without kids? The figure (below) shows an estimate for couple households without children where the household reference person is less than 55 years old. They really like the inner-west, dont they?

Back to my knitting

After a couple of forays into economic matters, I'll get back to the sort of stuff I normally do.

Below is a map (click the image for bigger version) showing the proportion of journey to work trips made by either walking or bicycle, at the time of the 2001 census. Not surprising to see the huge walk/bike mode share near the city.


Just FYI, the dark blue patch out in the north-west (near Richmond) is an army barracks, so lots of army dudes live near thhe barrackds and stroll in I guess. Same for in the south (near Holdsworthy). Other than that, its basically the CBD (incl. North Sydney) and the University of New South Wales, with Sydney's second CBD (Parramatta) also having a reasonable amount of walking and cycling.

More on Travel Behavior

Here are a few more graphs showing travel behaviour in Sydney.

The first shows average household VKT. This is estimated in the following way:
1) I estimate a simple household travel model that predicts household VKT based on household type, income, and location.
2) I take data from the 2006 census (which tells us where households are living) and estimate how far they will drive per annum

The model is pretty simple, but the trend is clear.


The above map is interesting, but there is a lot going on -- vehicle use changes with income and household type and location, and these all vary spatially.

Another way of looking at it is to think about the independent effect of location, so, what happens if you take the same household and move that household around. How does their vehicle use vary?

I cant do this exactly, for technical reasons I wont bore you with, but the map below is an attempt to get close. The actual method used to obtain this estimate is a little involved, and involves some guesswork, so dont take it as gospel, but it is at attempt to estimate the independent effect of location. In this instance, I take households from a wealthy area (but not stinking rich, i.e. not Vaucluse/Mosman) that has lots of 'couples with kids' households, and 'move' those households around (en-masse). The following map shows average household VKT when you hold income and household type constant, and change location only. I hope that explanation is comprehensible.


The last map (below) shows average commute time (by origin) in Sydney. This one is compiled as follows:

1) Take ABS journey to work data (They have origin -> destination trip counts by travelzone, which are the spatial regions shown here)

2) Take NSW Transport Data Centre zone->zone travel times (at peak hour).

3) From (1) and (2), calculate average commute per zone.

Again, nothing so surprising here.

Note that the commute time shown here is the average commute for all modes (i.e. car and walk and publilc transport and.. well... everything).