PerryGeo

In commemoration of the OGC approval of KML as an open standard to share geographic content over the web, I’d like to share our recent “Hike of Doom #2″ (kml provided by Mark Dotson).

The first weekend to hit 90 degrees, my friends and I travel inland to dive and swim in the Santa Ynez river. It is billed as a “30 minute” hike to our favorite watering hole. It becomes much more than that.

Of course the road leading up to the trailhead is closed due to construction so we have 3 miles of hiking on pavement just to get to the former trailhead- the Red Rocks parking lot.

Then the fun begins. A decent rainy season and some dam releases make for high flows and we’ve got half-dozen major river crossings to contend with. The recent fires added a good deal of organic matter to the river and the algae has bloomed accordingly. It is a wet, hot, rocky and slimy hike.

We make it to the swimming hole and enjoy the day. We dive, laugh, have a few beers.

The sun sets and the fun really begins.

Klaus, the Bavarian cyclist whom we’d met at the swimming hole, met up with us just after my girlfriend, Joselyne, sprained her ankle on a rock. Her ankle hadn’t started to swell yet but I could tell, drawing from my basketball injuries from the past, that she was not putting weight on it any time soon. We fashioned crutches from some driftwood. We met up with some turkey hunters (dressed in more camouflage more effective than most military uniforms) who helped us out by providing us some ankle wrap.

David and Andy began the trek back to the car to get help. The rest of us could either go back via the river bed , a rocky and treacherous endeavor given the setting sun, or head up to the main road and get some help. We decided on the main road and Shaun took off to alert the others to our plans. The main fire road was a trek in the opposite direction - longer, more elevation changes but smooth enough for a bike or truck and more accessible to vehicles.

I carried Jos, over my shoulder fireman-style and/or piggy-back, over the river crossings.

On the flats, Mark and I pushed Jos on Klaus’ bike.

We pushed on up the trail until we reached the main road. Klaus, after drinking the last of our beer, biked up to the dam keeper’s residence at Gibraltar Dam while Christina, Sarah, Mark, Jos and I continued up the trail. A half-hour later, Klaus and the dam keeper arrived in a pickup and drove the rest of us back to the Red Rock “parking lot”.

But the construction and rebar on the causeway meant there was no way to cross with a normal vehicle so we went by foot. Jos got back on Klaus’ bike and we pushed.

Luckily the slight downhill grade allowed her to glide back for a good portion, graciously sparing Mark and I from permanent back injury.

Meanwhile the away team had gotten some semblance of cellular reception and attempted to call the authorities. The goal was to get a ranger truck to drive out to get us or at least unlock the gate to meet us half way at the Red Rock parking lot. The authorities response was fantastic if not a bit overzealous. By the time we had gotten within a 1/4 mile of our car, we spotted helicopters. Then a firetruck. Then an ambulance. Joselyne was coasting by on Klaus’ bike and they didn’t even stop for her on the first pass! Apparently expecting to rescue a mangled body from the wilderness, the EMTs were somewhat disappointed at the less challenging situation they faced - a girl, coasting down the road on a bike with a sprained ankle.

We were back in the car, on the road before dark and got home in time for pizza.

So what did we learn from this? Well as a Boy Scout, I am ashamed to say I wasn’t prepared. A well prepped emergency kit would have helped a lot. At least we had an LED headlamp. Some rope would have gone a long way towards making a stretcher. An instant-ice-pack, ankle wrap and some ibuprofen would have been handy. We were wet and the mercury was falling quickly; some emergency shelter and clothing would have assuaged my concerns about the nighttime chill.

But this was offset by the generosity of the many people we met for the first time - The hunters who lent us their medical supplies, the dam keeper who got up from his Sunday dinner to make sure we got back safely, the EMTs who put tremendous resources into organizing a military-scale search party, Klauss who so generously stuck with us and shared with us his bike, his wisdom and his company. Without their help and our group of friends, the story might have a less happy ending.

Never underestimate the power of human kindness, generosity and cooperation! And never believe me when I say it’s a short hike.

I love python for its beautiful code and practicality. But it’s not going to win a pure speed race with most languages. Most people think of speed and ease-of-use as polar opposites - probably because they remember the pain of writing C. Cython tries to eliminate that duality and lets you have python syntax with C data types and functions - the best of both worlds. Keeping in mind that I’m by no means an expert at this, here are my notes based on my first real experiment with Cython:

EDIT: Based on some feedback I’ve received there seems to be some confusion - Cython is for generating C extensions to Python not standalone programs. The whole point is to speed up an existing python app one function at a time. No rewriting the whole application in C or Lisp. No writing C extensions by hand. Just an easy way to get C speed and C data types into your slow python functions.

So lets say we want to make this function faster. It is the “great circle” calculation, a quick spherical trig problem to calculate distance along the earth’s surface between two points:

p1.py

import math

def great_circle(lon1,lat1,lon2,lat2):
    radius = 3956 #miles
    x = math.pi/180.0

    a = (90.0-lat1)*(x)
    b = (90.0-lat2)*(x)
    theta = (lon2-lon1)*(x)
    c = math.acos((math.cos(a)*math.cos(b)) +
                  (math.sin(a)*math.sin(b)*math.cos(theta)))
    return radius*c

Lets try it out and time it over 1/2 million function calls:

import timeit  

lon1, lat1, lon2, lat2 = -72.345, 34.323, -61.823, 54.826
num = 500000

t = timeit.Timer("p1.great_circle(%f,%f,%f,%f)" % (lon1,lat1,lon2,lat2),
                       "import p1")
print "Pure python function", t.timeit(num), "sec"

About 2.2 seconds. Too slow!

Lets try a quick rewrite in Cython and see if that makes a difference:
c1.pyx

import math

def great_circle(float lon1,float lat1,float lon2,float lat2):
    cdef float radius = 3956.0
    cdef float pi = 3.14159265
    cdef float x = pi/180.0
    cdef float a,b,theta,c

    a = (90.0-lat1)*(x)
    b = (90.0-lat2)*(x)
    theta = (lon2-lon1)*(x)
    c = math.acos((math.cos(a)*math.cos(b)) + (math.sin(a)*math.sin(b)*math.cos(theta)))
    return radius*c

Notice that we still import math - cython lets you mix and match python and C data types to some extent. The conversion is handled automatically though not without cost. In this example all we’ve done is define a python function, declare its input parameters to be floats, and declare a static C float data type for all the variables. It still uses the python math module to do the calcs.

Now we need to convert this to C code and compile the python extension. The best way to do this is through a setup.py distutils script. But we’ll do it the manual way for now to see what’s happening:

# this will create a c1.c file - the C source code to build a python extension
cython c1.pyx

# Compile the object file
gcc -c -fPIC -I/usr/include/python2.5/ c1.c

# Link it into a shared library
gcc -shared c1.o -o c1.so

Now you should have a c1.so (or .dll) file which can be imported in python. Lets give it a run:

    t = timeit.Timer("c1.great_circle(%f,%f,%f,%f)" % (lon1,lat1,lon2,lat2),
                     "import c1")
    print "Cython function (still using python math)", t.timeit(num), "sec"

About 1.8 seconds. Not the kind of speedup we were hoping for but its a start. The bottleneck must be in the usage of the python math module. Lets use the C standard library trig functions instead:

c2.pyx

cdef extern from "math.h":
    float cosf(float theta)
    float sinf(float theta)
    float acosf(float theta)

def great_circle(float lon1,float lat1,float lon2,float lat2):
    cdef float radius = 3956.0
    cdef float pi = 3.14159265
    cdef float x = pi/180.0
    cdef float a,b,theta,c

    a = (90.0-lat1)*(x)
    b = (90.0-lat2)*(x)
    theta = (lon2-lon1)*(x)
    c = acosf((cosf(a)*cosf(b)) + (sinf(a)*sinf(b)*cosf(theta)))
    return radius*c

Instead of importing the math module, we use cdef extern which uses the C function declarations from the specified include header (in this case math.h from the C standard library). We’ve replaced the calls to some of the expensive python functions and are ready to build the new shared library and re-test:

    t = timeit.Timer("c2.great_circle(%f,%f,%f,%f)" % (lon1,lat1,lon2,lat2),
                     "import c2")
    print "Cython function (using trig function from math.h)", t.timeit(num), "sec"

Now that’s a bit more like it. 0.4 seconds - a 5x speed increase over the pure python function. What else can we do to speed things up? Well c2.great_circle() is still a python function which means that calling it incurs the overhead of the python API, constructing the argument tuple, etc. If we could write it as a pure C function, we might be able to speed things up a bit.

c3.pyx

cdef extern from "math.h":
    float cosf(float theta)
    float sinf(float theta)
    float acosf(float theta)

cdef float _great_circle(float lon1,float lat1,float lon2,float lat2):
    cdef float radius = 3956.0
    cdef float pi = 3.14159265
    cdef float x = pi/180.0
    cdef float a,b,theta,c

    a = (90.0-lat1)*(x)
    b = (90.0-lat2)*(x)
    theta = (lon2-lon1)*(x)
    c = acosf((cosf(a)*cosf(b)) + (sinf(a)*sinf(b)*cosf(theta)))
    return radius*c

def great_circle(float lon1,float lat1,float lon2,float lat2,int num):
    cdef int i
    cdef float x
    for i from 0 < = i < num:
        x = _great_circle(lon1,lat1,lon2,lat2)
    return x

Notice that we still have a python function wrapper (def) which takes an extra argument, num. The looping is done inside this function with for i from 0 < = i < num: instead of the more pythonic but slower for i in range(num):. The actual work is done in a C function (cdef) which returns float type. This runs in 0.2 seconds - a 10x speed boost over the original python function.

Just to confirm that we’re doing things optimally, lets write a little app in pure C and time it:

#include <math .h>
#include <stdio .h>
#define NUM 500000

float great_circle(float lon1, float lat1, float lon2, float lat2){
    float radius = 3956.0;
    float pi = 3.14159265;
    float x = pi/180.0;
    float a,b,theta,c;

    a = (90.0-lat1)*(x);
    b = (90.0-lat2)*(x);
    theta = (lon2-lon1)*(x);
    c = acos((cos(a)*cos(b)) + (sin(a)*sin(b)*cos(theta)));
    return radius*c;
}

int main() {
    int i;
    float x;
    for (i=0; i < = NUM; i++)
        x = great_circle(-72.345, 34.323, -61.823, 54.826);
    printf("%f", x);
}

Now compile it with gcc -lm -o ctest ctest.c and test it with time ./ctest… about 0.2 seconds as well. This gives me confidence that my Cython extension is at least as efficient as my C code (which probably isn’t saying much as my C skills are weak).

Some cases will be more or less optimal for cython depending on how much looping, number-crunching and python-function-calling are slowing you down. In some cases people have reported 100 to 1000x speed boosts. For other tasks it might not be so helpful. Before going crazy rewriting your python code in Cython, keep this in mind:

“We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil.” — Donald Knuth

In other words, write your program in python first and see if it works alright. Most of the time it will… some times it will bog down. Use a profiler to find the slow functions and re-implement them in cython and you should see a quick return on investment.

Links:
WorldMill - a python module by Sean Gillies which uses Cython to provide a fast, clean python interface to the libgdal library for handling vector geospatial data.

Writing Fast Pyrex code (Pyrex is the predecessor of Cython with similar goals and syntax)

Slashgeo pointed me to a very interesting set of projects - SpatiaLite and VirtualShape. They provide a spatial data engine for the sqlite database. Think of it as the PostGIS of SQLite. It looks like this extends sqlite’s spatial capabilities far beyond the sqlite OGR driver.

SpatiaLite provides many of the basic OGC Simple Features functions - transforming geometries between projections, spatial operations of bounding boxes, and some basic functions to disect, analyze and export geometries.

VirtualShape provides the really neat ability to access a shapefile using the SpatiaLite/SQlite interface without having to import a copy - it reads directly off the shapefile by exposing the shapefile and its attributes as a “virtual table”. I can think of a million uses for this. For example, lets say you have a shapefile of US counties and the number of voter in the 2004 election as an attribute in the dbf. You want to find the total voter count in each state:

$ ls -1 counties.*
counties.dbf
counties.prj
counties.shp
counties.shx
$ sqlite3 test.db
sqlite> .load 'SpatiaLite.so'
sqlite> .load 'VirtualShape.so'
sqlite> CREATE virtual table virtual_counties using VirtualShape(counties);
sqlite> select sum(voters) as total_voters, state_name
           from virtual_counties
           group by state_name
           order by total_voters desc;
9830550.0|California
7563055.0|Florida
7346779.0|Texas
...

Now this is fairly straightforward non-spatial SQL but the ability to run it against a shapfile without having to export to an intermediate data format is a very valuable tool.

Links:
When to use SQlite.
A video presentation by Richard Hipp (the author of sqlite).

What an odd meme .. I don’t know why but I expected some more interesting results. I guess the majority of the commands I use are pretty pedestrian.

history|awk '{a[$2]++ } END{for(i in a){print a[i] " " i}}'|sort -rn|head
163 vi
48 screen
29 python
28 ls
17 cp
17 cd
9 sqlite3
6 rm
5 sudo
4 htop

I don’t spend much time with web programming these days but I decided to give web.py (the minimalist python web framework) a shot and, while I was at it, try implementing a simple REST api.

First of all, web.py is truly everything it claims to be - small, light and easy to deploy behind lighttpd. It gives you a ton of flexibility to implement anything however you want - which is a plus or minus depending on how you look at it. I liked the inifinte flexibility but I can see alot of refactoring taking place and features needing to be implemented just to match the functionality built into a more structured framework like Django.

Back to the REST side of things. So I created a url-mapping to my “resources” or “nouns” and used the HTTP verbs (POST,GET,PUT,DELETE) to supply the interface. This was a joy to do in web.py which made it easy.

urls = ("/thing/(\d+)", "thing")
...
class thing:
    def GET(self, thingid):
        # select query and render to template
        ....
    def POST(self, thingid):
        # insert query and redirect to /thing/thingid
        ....
    def DELETE(self, thingid):
        # delete query
        ....
    def PUT(self, thingid):
        # use cgi args to run update query on specified thing
        ....

The hard part came when I realized that HTML forms do not implement DELETE or PUT methods! 2 of the 4 cornerstone HTTP verbs are not implemented in HTML forms?

Surely this can be accomplished with a top-notch AJAX library. I tried Prototype.js and it appears that the PUT and DELETE methods are simply tunneled over POST with an extra arg attached and the server side has to handle it accordingly. So I ended up just using a straight XMLHttpRequest which works but has it’s own problems.

How are you supposed to call PUT or DELETE through a web page? Is XMLHttpRequest the only way? What about browsers without javascript?

There are two new books coming out this summer which fill a valuable niche in the open-source GIS bookshelf:

• Applied Spatial Data Analysis with R
• Desktop GIS: Mapping the Planet with Open Source

These are both written by some of the top developers within their respective topics and I’m really looking forward to reading them.

The X61s is one bad-ass machine. Besides the great performance, battery life and solid engineering, there are other hidden gems. Like the tilt sensors that were designed to protect the hard drive in case of a drop can also be used to detect the laptops motion under more normal circumstances.

There are some interesting applications that use some simple statistics to determine when the machine is “tapped” or julted to left or right. You can then assign actions to unique combinations of taps.

These applications all use the sysfs interface to the sensors ( cat /sys/bus/platform/devices/hdaps/position will show your position in the x and y axis). But the sensors also provide a joystick interface that allow you to tilt the laptop along the two horizontal axes to control any number of applications. Including Google Earth.

1. Install tp_smapi
2. Test the sensors by running hdaps-gl , a simple OpenGL app showing the real-time tilt of your thinkpad.
3. Run jscal to calibrate the joystick. You’ll need to install the “joystick” package for this. The command is:
   jscal -c /dev/input/js0
   After which you should keep your laptop level for a few seconds. Then, when prompted, tilt left, center, right, back (towards you), center, then forward.
4. Now fire up Google Earth. Open the Options menu, go to Navigation and select Enable Contoller.



5. You should now be able to zoom around by tilting the laptop. The keyboard shortcuts really help when you’re in this mode (Ctl-Up/Down to zoom, Shift-Up/Down to tilt, Shift-Left/Right to pivot).

There’s also a neat Perl-script technique to control a web-based google map which has some cool potential for an openlayers based system.

Since most Apple laptops have a similar sensor, you should be able to get the same thing going on your Macbook. Try it out..its alot more fun that using the mouse!

My HP laptop was nearing 5 years old. It had held up extremely well but most modern software taxed it to the absolute limits (just having firefox open with a flash ad in one tab was enough to send the system load through the roof). So I decided to try something new.

I was looking for something in the ultra-portable range. I tried out the OLPC and looked seriously at the Asus eeepc for a while. But they were far too difficult for me to type on. Ergonomics were extremely important and the only ultraportable that consistently rated high in that department was the IBM/Lenovo thinkpads. The X61s was appealing with its low voltage core2 duo and 2GB of RAM. All that in a small package about 3 lbs and about an inch thick.

So the X61s arrived and I figured I’d give it a try with the “stock” software. It was my first experience with Vista and I gave it my best shot. After about 1/2 hour of excessive clicking, boggy performance and pop-up windows, I shrunk the ntfs partition and installed Ubuntu Hardy Heron Alpha 4.

Sound, wireless with WPA, Compiz with 3D; the major things that normally plague a linux laptop install worked right out of the box. On the other hand, I’m running into a few bugs in nautilus (this is is alpha software after all), I can’t get bluetooth working, suspending to ram works but is a little buggy (have to restart some services manually) and I had to edit a few config files and compile a kernel module to utilize all the bells and whistles provided by the hardware. But it is still more fun than using Vista.

One thing that really shines on this machine is the battery. I got the 8-cell extended life battery and used some powertop tweaks cut my power consumption and was able to get the wattage down in the 10 to 15 watt range depending on usage patterns. No wonder it is energy star compliant! With that kind of wattage and battery capacity, I’m easily getting about 6 to 7 full hours of battery life.

Some tips if you’re setting up Linux on your X61s:

We did it!

As some of you may know in 2005 through 2006, I was part of a research team[1] , led by Ben Halpern at NCEAS, developing a global model of human impacts on the marine ecosystem. We created or compiled 17 high-resolution global datasets of human-induced threats (land-based pollutants, fishing, shipping, climate change, etc.) and 20 ocean habitat datasets. These were combined to create an impact index which models the cumulative level of human-induced stress on our oceans.

The results were published today in Science magazine and presented yesterday at the AAAS Annual Meeting. To summarize, we found that the entire ocean is affected and 40% is heavily impacted. It is not all bad news as there are many areas of relatively low impact which could provide examples for ecosystem restoration and opportunities for conservation. The global map is the first of its kind and will help clarify and quantify our cumulative impacts on the ocean and allow us to focus efforts geographically. The model is not perfect and can’t really be used to make decisions at a very localized scale but, given the available globally-consistent, reasonably-high-resolution data for all the various ocean threats and habitats, this is the best effort to date. The model itself is relatively simple with a very clear methodology which will allow scientists to tweak the parameters and add better data as it becomes available. For those of you interested in the GIS modeling end, NCEAS has a great summary of the data used in the model. Most of the data are available as raster data products or KML.

The media has picked up on the story with NPR, MSNBC, The Washington Post, USA Today and National Geographic covering it (to name a few). I especially recommend the NPR site as it has a great animation and an audio segment.

So congratulations to everyone who made this happen!

[1] Benjamin S. Halpern, Shaun Walbridge, Kimberly A. Selkoe, Carrie V. Kappel, Fiorenza Micheli, Caterina D’Agrosa, John F. Bruno, Kenneth S. Casey, Colin Ebert, Helen E. Fox, Rod Fujita, Dennis Heinemann, Hunter S. Lenihan, Elizabeth M.P. Madin, Matthew T. Perry, Elizabeth R. Selig, Mark Spalding, Robert Steneck, Reg Watson (2008). A global map of human impact on marine ecosystems. Science, vol. 319

EDIT:

Some additional articles:

• New York Times
• BBC Video on YouTube

Sadly, a lot of GIS folks have never come into contact with a command line interface (CLI) . I’ve met even experienced computer users who, when faced with a command-line prompt, experience some autonomous nervous system lock up that causes their eyes to glaze over and prevents any knowledge from entering their brain from that moment forward. The all-Windows, all-GUI mentality of the current GIS market leaders just doesn’t expose you to it (if you remember working with coverages at the ESRI Arc/Info command line, you official qualify as an “old-timer”). And the DOS command line is virtually invisible to XP and vista users. Linux users are more CLI aware but this is even becoming less important as distros such as ubuntu GUI-ify everything.

So why the fear of the command line? Why is it assumed to be more “complicated” than a graphical user interface (GUI)? I have found that, in some cases, the opposite is true … there is something reassuringly simple about typing something and getting a response back. It feels like you are in direct control of the computer. Which, indeed, you always are. The computers always do exactly what you tell them, whether you are in a GUI or a CLI. But GUIs attempt to abstract away the details so that you don’t need to know exactly what you’re telling the computer to do. This nice fluffy feeling comes at the cost of many important factors. Consider:

• Automation: If you had, for instance, monitoring data coming in in a hourly basis and needed to process the data, would you want to be on call 24 hours a day to click a few buttons. Of course not. Write a command that performs the job and schedule it to execute at some regular interval. (I wonder if those guys on LOST ever thought to just set up a cron job to enter the numbers in the hatch?)
• Repeatability Whenever I show someone a CLI-based method for solving their problem, they almost immediately say (or at least imply) that the typing is too much trouble. Consider this command to convert a .tif image to ERDAS .img (HFA) format:

  cd /data/images
  gdal_translate -of HFA aerial.tif aerial.tif.img

  You might ask, “Why not just use a GUI, click a button or two, and get your output”. Sure. Now do that for 2,000 tif images. With a CLI you only have to type a few extra lines.

  cd /data/images
  for i in *.tif; do
    gdal_translate -of HFA $i $i.img;
  done

• Documentability: There is nothing more important to a GIS Analyst than documenting his/her work! We live by metadata and methods write-ups. Now picture an intense 5 hour work session … everything needed to get out by 2pm. You’re done and now it’s time to document your procedure and methods. With the CLI, you copy and paste your commands from the terminal or simply look at your command history which will show exactly what you did and how. You can store this in a text file and come back to it months later and be able to re-run the procedure.

  With the GUI, you have to remember and describe every click, every sub-menu, every option, every action taken to arrive at the answer. Often this requires verbose description, screenshots, etc. None of which is recorded in any history file of course. And of course, when the client inevitably comes back the next day with modifications, none of it is repeatable in any automated fashion with a GUI.

• Accessibility: It’s just plain text with a CLI. You can print it out and study it on the bus. You can email the whole process to co-workers. You can use a concurrent versioning system to keep track of changes to scripted procedures. You can transfer massive amounts of knowledge without having to sit down and go through everything step-by-step, click-by-click in a visual interface.
• Accuracy: Far too often, GUI designers make over-reaching assumptions about how things should work. The idea is often that the user should not need to know anything more than the absolute minumum. To use a car analogy, the driver turns the key, presses the pedal and steers but does not need to know what goes on under the hood. This works most of the time. But the law of leaky abstractions usually takes hold and something inevitably breaks or performs differently than expected. Since the CLI does not hold your hand (it executes the exact command you give it) it more accurately mimics the actual physical interaction with the computer and is much more useful in debugging and investigating complex problems.

So basically, don’t make the mistake of thinking that a pretty window will always contain the magic button to get the job done. In many cases, a command line is much more efficient, even essential. If you don’t know how to effectively work in a command-line environment, do yourself a huge favor and learn.

Oh and I’d be remiss if I didn’t mention Neal Stephenson’s book on the subject … a bit technically outdated but a great quick read on why command lines are still very relevant in the face of increasingly sophisticated graphical interfaces.