Has Automation affected YOUR job?

Can't resist posting this (from Modern Mechanix).

I'm definitely enjoying not having grimy hands, but what happened to the shorter work week and more leisure time?

World Population Cartograms

The Daily Green has some cartograms showing the change in distribution of world population since 1900.

But there's another aspect to this that the original images don't show - the increase in the absolute number of people in the world. So here's the images with that factor applied:

1900 - 1.65 B

2007 - 6.64 B

2050 - ~ 9 B

A good kick in the pants for the GeoWeb?

Reading the blog posts that are starting to come out about the ESRI Dev Summit, I'm struck by a few things:

• REST! KML! JSON! AJAX! Dojo! I haven't seen this many tech buzzwords out of ESRI since ArcGIS on Windows was released with COM! VBA! DLLs! Hopefully some of these have more legs
• Clearly ESRI is embracing the GeoWeb concept with a vengeance. It will be very interesting to see how these new modes of access to GIS functionality plays out. To what extent will the marriage of web services and GIS processing be an effective model?
• They also seem to be accepting Google Maps/Earth and MS VE as valid spatial delivery platforms. It sounds like there is some pretty full-featured capabilities to deliver data to these platforms. It's going to be interesting to see the uptake on these capabilities, and whether ESRI's heart is really in making these platforms perform to their full capability. No doubt there will be some fascinating business implications as well that get played out over time...
• I'll be interested to see whether ESRI's "legitimizing" spatial-via-web-service will have any effect on the OGC W*S world. It seems to me that while OGC was early out of the gate with a web service suite, they seem to be wallowing in the doldrums as far as making the specs effective for real-world use. The OGC W*S suite was quickly adopted by the OSS world, and for good reason - open source tends to be very faithful to open standards, since a design goal is usually to have a high degree of interoperability. Not to mention that it's easier to code to a standard which has already done a lot of the hard thinking. But - the cool kids are losing interest in the crusty old W*S interfaces, since they're not keeping up with the rapid emergence of exciting new web paradigms.
• Ultimately ESRI's trailblazing will be a good thing for OSS technologies. Nothing like having a working system to inspect, copy, and improve. And the more open, standard technologies that system uses, the easier it is to steal - um, learn from. Also, as the GIS stack gets more open interface technologies accepted, it becomes easier to plug in heterogeneous components into that stack.

Now to start looking into that ArcGIS Server REST API...

Tired of those boring old conformal map projections? Make your own!

Flex Projector is an interesting tool that allows you to make a custom map projection by adjusting the location and shape of meridians and parallels.

Here's a nice projection I defined to maximise the area of tropical vacation spots 8^)


The downside of this app is that as far as I can tell, there's no way to export a projection once it's defined. Instead, you have to use the Flex Projector application to reproject datasets. So it doesn't look like you can make your own custom projection for use in say, PostGIS or OGR.

I guess it might be asking a little much for the app to spit out some custom Proj4 C code with an associated EPSG code and CRS WKT... But hey, it's an open source app, so come on someone - how about it?

Quote of the (Pi) Day

Sir, I send a rhyme excelling
In sacred truth and rigid spelling
Numerical sprites elucidate
For me the lexicon's full weight

Maybe not the greatest lyric poem - but list the number of letters in each word. (If you prefer a shortcut this site has them already listed - along with 9, 979 more)

And notice the value of the date as MM.DD...

Branch-and-Bound algorithms for Nearest Neighbour queries

This paper by Roussopoulos et al is a fine expositions of how to use a Branch-and-Bound algorithm in conjunction with an R-tree index to efficiently perform Nearest-Neighbour queries on a spatial database.


Nearest-neighbour is a pretty standard query for spatial database. Oracle Spatial offers this capability natively. It would be great if PostGIS did too. If the GIST index API offers appropriate access methods, it seems like it might not be too hard to implement the Roussopoulos algorithm. (How about it, Paul, now that you're dedicating your life to becoming a PostGIS coding wizard...)

I'm also thinking that this approach might produce an efficient algorithm for computing distance between Geometrys for JTS. It's always bugged me that the JTS distance algorithm is just plain ol' O(N^2) brute-force. It seems like there has to be a better way, but as usual there seems to be precious little prior art out there in web-land. This should also work well in the "Prepared" paradigm, which means that it will provide an efficient implementation for PostGIS as well. Soon, hopefully...

The Four Programmers

Here's an amusing adaptation of Monty Python's Four Yorkshiremen sketch...

I feel increasingly like this myself (more so now that I'm in a company where I am the oldest employee)...

When I started coding we had to type out our programs on punchcards. If you made one typo you had to retype the entire card, so to avoid mistakes we wrote out our code on paper pads with 80-column grids. You had to stand in line to submit your card stack to the guy running the card reader, and then pick up your output once it had been separated and filed by some other guy. And if you had a syntax error - back to the cardpunch to do it all over again!

We thought we'd died and gone to heaven when we got online accounts and were allowed to use a hardcopy DecWriter terminal. It even had an APL character set (don't forget to change the charset when you switched back to FORTRAN). This was better - we could "erase" (actually just strikeout) and change mistakes on the line. But forget about printing out your entire program - at 300 baud this could take many minutes, which was all chewing up your connect time allocation.

But we had it lucky.. I worked with a guy who started out coding accounting routines for a magnetic drum-based system in the early 60's. A big chunk of their time was spent reordering the individual instructions around the drum to reduce read latency. Fun stuff!

And he was lucky compared to Don Booth, who was at the oceanographic research centre where I worked during my undergrad days. He was one of the pioneers of computing in Britain. All they had for storage was a 1000-word mercury delay tube... and no doubt they spent their coding time plugging wires and changing tubes.

But you try and tell kids with their touch-sensitive palmtops that these days, and they won't believe you...

The End of an Architectural Era?

Stonebraker et al have an interesting paper pointing out the antiquity and consequent limitations of the classic relational database architecture in today's world of massive disk/cycles/core.

If they are correct (and in spite of the recent MapReduce blunder Stonebraker has made a lot of great calls in the DB world), the world of data management is going to get awfully interesting in the coming years. The DBA's & DA's of this world have been living a relatively comfortable existence compared to those who are wandering across the stormy badlands of the middle tier. But Stonebraker postulates at least 5 radically variant database architectures to address specific use cases of data management. This would seem to lead to a much more complex world for data architects. But maybe a windfall for the consultants who find their niche?

Quote of the Day

Java - learn once, write anywhere

"From my personal experience, that is quite true. I started programming in Java in 1997, and since then I have built web applications, database applications, mobile applications for Palm and Pocket PC, an application server for CTI applications, CORBA-based distributed applications, network management applications, server-side business components for application servers, fat-client business applications, an IDE, a Rich Client Platform, a modeling tool among other kinds of software. All that across different operating systems and using similar products by different vendors.

Of course I had to learn something new for each kind of "application style", but the programming language and core class library were always there for me. That is a huge thing. Before Java, it was really hard to move from one domain to another, or even move from vendor to vendor. Sun has my eternal gratitude for bringing some sanity to the software development industry."

- Rafael Chavez

The Paradoxically Poor Power of Visual Programming

I've just been reading about TextUML, which is a language and IDE for textually defining UML diagrams by a start-up company called abstratttechnologies.

This seems like both a stroke of brilliance and painfully obvious. It makes me wonder why UML isn't already defined in terms of a serializable language. The big advantage of this I think would be the ability to have a clearly defined semantics. In the past when reading about UML I've been quite frustrated by an apparent lack of rigour in the definition of the various UML concepts and diagrams. It's always seemed to me that a metalanguage like UML is more in need of fully-defined semantics than application languages. Most (all?) application languages after all are executable, which gives provides a final authority for the meaning of a program written in them - run it and see what happens. AFAIK UML isn't really (or not conveniently) executable in the same way.

Which brings me to the paradox promised in the title of this post. One of the claims for TextUML is that it is easier to write and express complex models in than a purely diagram-based UML tool. I can easily believe this - developing a complex model by diagrams seems cumbersome and limited compared to writing a textual description. And I think this generalizes to the concept of visual programming in general. We've had easy access to powerful graphics capabilities for over 20 years now, culminating in today's basically unlimited in graphical capabilities. In spite of this, text-base programming still rules the roost when it comes to expressiveness, flexibility and speed.

So the paradox is: why is a 1-dimensional representation more powerful than a 2-dimensional representation? To put this another way, why does a 1D representation seem less "flat" than a 2D representation?

For further study: would the same apply to a 3D representation? Or would a move to 3D jump out of some sort of planar well, and end up being more expressive than 1D?

By the way, abstratt is right here in l'il ol' Victoria, and the lead designer (and CEO) Rafael Chavez is giving a talk at VicJUG in March. I'm definitely going to attend - I suspect he may have an interesting take on this issue.

New Buffer Styles in JTS 1.9

JTS allows specifying styling parameters to control the appearance of generated buffer polygons.

End Cap styles have been available for a while. End caps can be Round (the default), Flat, or Square. For Round caps, the Quadrant Segments parameter controls the degree of approximation to a true circular arc.

End Cap Styles: Round, Flat, Square

JTS 1.9 adds the ability to specify the Join Style. The Join style controls how the offset curves for consecutive line segments are joined together. The Join type can be Round (the default), Mitre, or Bevel. If Mitre is specified, a further Mitre Limit parameter is provided to control how long mitred joins can be before they are beveled off. The maximum length of a mitred corner is equal to the Mitre Limit multiplied by the buffer distance.

Join Styles: Mitre (limit = 10), Mitre (limit = 1), Bevel

The JTS TestBuilder has been enhanced to allow specifying all of these parameters, so it can be used to experiment with the visual effects of the various styles.

Hey What's That - Terrain Panoramas Mashup

Here's a nifty mashup using Google Maps. It computes a skyline terrain panorama from a given location.



The computed skyline isn't bad, although a bit optimistic (I'm pretty sure you can't see the Olympics from downtown Vancouver!). I think this particular query is ignoring earth curvature (although apparently you can choose to take it into account.)

The site uses it's own SRTM dataset. The gazetteer is impressive - not sure how they get that information.

Cross-Border issues with Google Maps terrain

So GM now has terrain - yay! But they seem to have short-changed half the continent. Check this screenshot from the border just south of Chilliwack. Notice the American Border Peak - but where's the Canadian one? And what's with that crazy-looking border swathe?


I assure US viewers that Canada isn't normally this blurry. And no, it isn't proof that we really do have a year-round 4 ft snowpack!

Ok, I know - never satisfied. Really, this is a very nice feature. I just hope Google gets their north-of-the-border vision lasered real soon...

Enumerating permutations using inversion vectors and factoradic numbering

Permutations are a fundamental concept in discrete mathematics, and crop up in numerous places in computer science. The set of all possible permutations of n items is usually denoted S_n (since the set forms a group which is isomorphic to the Symmetric group of order n). The number of permutations in this set is n!.

This apparently simple concept poses some unexpected challenges. For instance, the following useful operations are messy to compute via the brute-force way of enumerating successive permutations:

• List all the permutations in S_n (in lexicographic order)
• Find the i'th permutation in S_n

This would be straightforward if there was a way of effectively enumerating all permutations; i.e. a mapping between S_n and the integers {1, 2, ..., n!}.

It turns out that this can be done using a pair of ingenious concepts known as the inversion vector and factoradic numbering. These can be used to provide an enumeration of S_n, conveniently in lexicographic order. (Inversion vectors are also known as inversion tables and Lehmer codes - the term vector seems to me to be most appropriate for the nature of the structure.)

A permutation can be written as a sequence of numbers

P = (p₁ p₂ ... p_n) where p_i is in {1, 2, ..., n}

An inversion is a pair (p_i, p_j) where i < j but p_i > p_j. We define d_i to be the number of inversions in the sequence where the first element is p_i:

d_i = | { p_j : i < j , p_i > p_j } |

Then the inversion vector for the permutation is the sequence

D = (d₁, d₂, ..., d_n) (d_i is in the range {0, 1, ..., n-i})

It is fairly straightforward to recover the original permutation from the inversion vector. This paper gives a definition of the forward and reverse mapping algorithms.

There are n! possible sequences D. They are in one-to-one relationship with the set of permutations S_n, and their natural order gives a lexicographic ordering of S_n. So we're halfway to the enumeration that we are looking for.

Inversion vectors can in turn be mapped to the integers in {1, ..., n!} by using the concept of factoradic numbering. A factoradic number is a number made up of a sequence of digits where the radix for digit d_i is (n-i)!. This Wikipedia article explains the concept in detail.

Here's an example:

P = (4 6 2 5 3 1 8 7) is a permutation in S₈.

Its inversion vector is L(P) = [3,4,1,2,1,0,1,0].

The factoradic representation of this is

3×7! + 4×6! + 1×5! + 2×4! + 1×3! + 0×2! + 1×1! + 0×0! = 18175

So inversion vectors and factoradic numbering provide an elegant way of enumerating permutations. But what does this have to do with geospatial processing? That will be a subject for another post...

If programming languages attended computer conferences

If you're a language geek, this post by David Rupp is hilarious!

Fast polygon merging in JTS using Cascaded Union

A common requirement in spatial processing is to union a set of polygons together. Depending on the use case, the polygons may either be a coverage or be overlapping. Up until recently, there were two techniques for accomplishing this in JTS:

• Iterated Union: Iterate over the polygons and union them one-by-one into a result polygon
• Buffer Union: Collect the polygons into a GeometryCollection and run buffer(0) on the collection

The first one is straightforward, but quite inefficient for large input sets. The second technique is more efficient in many cases, but still tends to fail for large or complex data.

Since this is such a common operation, I've been keen to find a faster and more robust approach. An alternative strategy is to use a concept I call Cascaded Union. The idea is to union small subsets of the input together, then union groups of the resulting polygons, and so on recursively until the final union of all areas is computed. This can be thought of as a post-order traversal of a tree, where the union is performed at each interior node. If the tree structure is determined by the spatial proximity of the input polygons, and there is overlap or adjacency in the input dataset, this algorithm can be quite efficient. This is because union operations higher in the tree are faster, since linework is "merged away" from the lower levels.

This picture shows the algorithm in action:


Handily, the spatial indexes already in JTS can be used to determine an appropriate tree structure. Currently the Sort-Tile-Recursive packed R-Tree index is used - this seems to produce an effective tree (although I suspect that the algorithm is not very senstive to the precise spatial tree used)

It turns out that Cascaded Union can be much more efficient than Iterated Union. On a test dataset with 30,000 small polygons with a high degree of overlap, the performance results were:

Cascaded Union: 20 sec
Iterated Union: 3 hours 40 min !

Not bad - well worth the effort of coding... Also, it's more robust than Buffer Union, and often much faster as well.

This will be provided via the Geometry.union() (Unary union) method in JTS 1.9.

Implementations of Polygon Overlay algorithms

Occasionally people ask about the algorithms used in JTS's overlay operations. I don't have a good answer for this, since I basically made it up (in an informed kind of way), and haven't had time to document it to the level I'd like (thank goodness for self-documenting code... 8^)

I recently came across this web site from back in 1998, which is happily still around. It has a good list of implementations of what it calls Polygon Boolean operations. (This is the same thing as the JTS overlay operations - I chose to use a different term to avoid confusion with the spatial predicates). Perhaps some of these sites have done a better job of documentation! There's also some useful references to academic papers covering aspects of the issue.

The bible of Spatial Indexing

The Bhagavad Gita of Multidimensional Access Methods...
The Koran of Quadtrees...

One of the people I met at ACM GIS 2007 was the conference chair, Hanan Samet. He's published a series of highly useful books on spatial indexing methods, which are certainly the most comprehensive and detailed in this field. His latest work is Foundations of Multidimensional and Metric Data Structures and is surely his magnum opus. It could also be called Everything you wanted to know about Spatial Indexes (but were afraid to code). Recommended reading for spatial geeks, especially if you're looking for some reading to tide you over during your next year off.

His first book on the subject, The Design and Analysis of Spatial Data Structures, is also worth looking at, especially since it's a bit more bite-sized. However, it's out-of-print, and expensive when you do find a copy (local plug: I found mine through ABE) If you really can't find a copy you might try contacting Hanan directly.

Hanan has a clear preference for quadtree structures over R-trees, but both approaches are covered well (including pseudo-code - yay!)

Conference Review: ACM GIS 2007

I just got back from ACM GIS 2007, over in Seattle. This is the first time I've attended this conference - or even known about it! It was a good year to go, since there was a lot of energy in the air due to this being the first year that the conference has been held as an independent event in its own right (which may be one reason I haven't heard of it before). One reason for doing this is to catch the building geospatial wave - more on this below.

The conference was very well attended (about 200 people). The bulk of attendees seemed to be from US universities - but hailed from all over the globe. Microsoft was there in force, as you might expect, and Google was very prominent too (a running joke at the conference was that every talk was connected to or given by someone tied to MS or Google). The omnipresent Michael Jones gave a good motivational talk about some of the technologies which are challenges for Google and would be fruitful research areas. Clearly the big guns of the GeoWeb are investing heavily in geospatial research, and this in turn is injecting a lot of cash and relevancy into this field.

Some of the presentations which stood out for me were:

• Using Constraint Programming to conflate address data to urban imagery. Constraint programming is a well established "AI" technique with a lot of mature effective algorithms, and it seems like a promising approach for other kinds of conflation problems as well (such as road network conflation)
• A guy from Google talked about conflating vector road networks to ortho imagery (surprise, surprise - this is a big use case for them). His approach seem to take place mostly in the raster domain (using edge detection and signal isolation), which I found interesting
• A talk by Andrew Danner about building massive TINs and extracting drainage patterns from them. Having just tackled a very similar problem (in both domain and data size) I was interested to see his approach - which was a fair bit more clever. His group has a novel external memory algorithm for building TINs on huge datasets. This seems like an important focus for research, since data gathering bandwidth is outpacing RAM growth.
  
• A couple of talks on approaches for compressing surface models and extracting varying levels of detail. Nice because they used some elegant mathematical concepts such as Morse-Smale codes.

There were lots of other talks - and a HUGE bonus was that the conference proceedings were provided to you at registration time (not 2 months later). So you could read ahead to get a good idea of what was coming up - or use downtime to get a more in-depth understanding of a previous talk.

Speaking of which, this is the first conference I've been to which had a single track of presentations. This made for pretty full days, and was probably only barely possible due to the size of this conference, but it certainly made for a more satisfying, less stressful, and more inclusive experience. No need to agonize on which one of 5 presentations might be the most interesting and then rush madly around repositioning during the 5-min break between talks. And it gave you time to catch up on your reading if the subject wasn't of total interest.

Next year the conference is in Irvine, CA. The organizers are keen to get more industry participation, so don't be shy - check it out.

Workaround to Google 4 Black Screen Of Death problem

Ever since I upgraded to Google Earth 4, I've had the extremely annoying problem that tilting to an oblique view produced a completely black screen. The only way to get the screen view back was to restart GE. Time to dust off that ol' BSOD acronym...

I just found a thread which provides a clue to why the problem occurs, and gives a workaround (which doesn't 100% remove the annoyance, since it only works for that one session).

Apparently the problem is caused by the use of a dual core AMD Athlon processor. The solution is to set the processor affinity to use only a single core for the GE process. (TaskManager, Processes Tab, rightclick process 'googleearth.exe', 'Set Affinity', select 1 processor only). I haven't yet found a way to persist this setting - anyone have any ideas?


Update: After a few milliseconds more browsing, I found that you can use the imagecfg.exe utility to permanently set the processor affinity for a load image. Running the following command in the Google Earth directory does the trick:

imagecfg -a 0x1 googleearth.exe