The Roshambo-bots are coming!

The Japanese have figured out a way to let me avoid taking the garbage out on rainy winter evenings - a Rock-Paper-Scissors robot that always wins!

It actually cheats, though - it reads your mind to determine what move you're going to throw.

Now, enough playing with the easy stuff - get back to work making a Go program!  Or at least Rock-Paper-Scissors-Lizard-Spock...

A scientific basis for Open Source Software

Stefan Steineger of the OpenJUMP project pointed out this great paper in Nature on The case for open compute programs.  The paper raises the argument for open source software to a higher plane, that of being a necessary component of scientific proof.  It points out that the increasing use of computational science as a basis for scientific discovery implies that open source must become a standard requirement for documentation.  Apparently some journals such as Science already require source code to be supplied along with submissions of articles.  Amongst other advantages, access to source code is an essential element of peer review.

An interesting example they mention is the infamous HadCRUT and CRUTEM3 meteorological datasets.  One of the (few) salient criticisms levelled at this information during Climategate was the inability to  reproduce the results by re-running the software. (Mind you, the software was probably a pile of crufty old Fortran programs mashed up by Perl scripts, so maybe it's just as well).

I'm looking forward to seeing JTS get cited in academic papers (actually, it already has been).  Maybe I even have a finite Erdos number!

It's maybe too much to ask that mere scientists be coding hipsters, but I noticed that SourceForge is presented as the leading example of collaborative software development.  Someone should introduce them to GitHub - which truly walks the talk.  Researchers in bioinformatics should be especially appreciative of the sweeping effect of recombinant software development it enables.

Why we were stronger in the 80's

I love this picture. Nowadays this functionality weighs only a few ounces, and we carry it in a pocket.

Grandpa's Googler

A peerless pastiche by SWWTMTOTH...  Try it, it works!  Go ahead, make it your homepage!

Hands up all those who remember modem whistles.  Takes me back...    For maximum verisimilitude, every other time you hit this site it should give you a busy signal.

Little known fact - DSL and cable modems actually make the same sounds, they're just pitched too high to hear  8^

FOSS4G-NA 2012 review

Last week I was at FOSS4G-NA in Washington DC.  It was my first time in DC, and my first FOSS4G as a member of the OpenGeo team. Both were very pleasant experiences!

One of the keynotes was by Josh Berkus of PostgresQL talking about "firehose databases".  FOSS4G is always a firehose conference, and this one was no exception.  I missed many more talks than I wanted to take in, for various reasons - the 3 simultaneous tracks, my own commitments to talks and OpenGeo work, and often because it's just such a good chance to talk to great people from all around the FOSS4G world.

Here are some of the my conference highlights:

• another inspiring keynote from Paul Ramsey, on the origins and future of open source
• the Ignite Spatial talks, which were fast, fascinating and funny
• a panel discussion which hit on the need for a replacement for the shapefile (PLEASE!)
• hearing Josh Berkus talk about the impressive new features coming in PostgresQL 9.2 (JSON, index-only scans, SP_GIST, to name a few)
• Andrew Turner's fastpacedpresentationonGeoIQ
• meeting Kevin Web and hearing more about OpenTripPlanner
• discussing skeletonization, long-distance hiking, and sewage treatment with Stephen Mather 
• talking to Josh Marcus about GeoTrellis, an engine for scalable geoprocessing
• talks on using NoSQL DBs such as MongoDB and HBase for storing and querying spatial data
• a discussion with Javier  de la Torre of CartoDB about the need for a spatially-sensitive sampling technique (a la Oracle's SAMPLE clause) in PostGIS.  We ended up discussing this with Josh Berkus (since he was sitting at the next table!) and apparently there is some work starting on this in PostgreSQL.  Hopefully this will extend to spatial datasets as well.
• and of course, meeting lots of people using JTS!

Here's a few other reviews from around the web:

• Directions Mag takeaways from the ever-present Adena Schutzberg
• smather on a Real American Hero
• Bill Dollins also had the firehose feeling
• Josh Berkus liked his sip from the spatial firehose

Battle of the SSLs - Round 2

In a comment to Battle of the Spatial Scripting Languages, Sean suggested that a slightly more in-depth test of language expressiveness might be this task, involving a sequence of CRS transformation and geometry cleaning.  Here's the JEQL equivalent:

ShapefileReader t file: "test_uk.shp";

trans = select * except GEOMETRY, 

               CRS.transform(GEOMETRY, "epsg:4326", "epsg:27700") geom from t;

clean = select * except geom, isValid ? 1 : 0 isValid, cleanGeom
    with { 
        isValid = Geom.isValid(geom);
        cleanGeom = isValid ? geom : geom; //Geom.buffer(geom, 0);
    } from trans;

ShapefileWriter clean file: "test_uk_valid.shp";

The code is spread out a bit for clarity.  It could have been done in single SELECT statement, apart from the isValid flag added for reporting purposes.  Note the use of the JEQL extension with clause, which improves code clarity, and also allows a more imperative style of coding which can sometimes be preferable.

Fiona gets points for being able to define any function directly in the language, but there's a price to be paid in code complexity.

Here's the result (seen in the new JEQL Workbench Geometry Viewer):

The TI-58/59 programmable calculators

Mobile computing is nothing new!  My first programmable device was a TI-58 calculator:

I spend many happy hours punching in the code for Lunar Lander (and less time actually playing it).  Wish I still had my coding sheets for the various apps I dreamed up.  (For some reason, the convention then was to write the key codes vertically, rather than horizontally.)

TI had their own app store, selling Solid State Software^TM on plug-in ROM modules.  Of course I got the Leisure one, which had a few lame games on it.  Then as now, though, I had a lot more fun coding up new apps than playing games...

The calculator cognoscenti will notice the image is actually the TI-59, with the built-in magnetic card reader for offline permanent storage.  Had one of those too, once I got tired of eternally punching in programs. The mag card reader was pretty finicky, but it sure was cool filling up that little vinyl case with all that software that never got used...

A couple of years later HP released the iPhone of the calculator era - the HP-41C, which had an alphanumeric display!  HP was just as opinionated as Apple, in their own way - if you wanted one of their gorgeous units, you had to grok RPN notation. Kinda reminds me of my current effort of learning git.  I never did get one of those units, because by that time it was a lot more rewarding to punch code into a real microcomputer....

Reminiscing about the PC Pre-Cambrian explosion

With the mobile meteorite looming ever larger over the landscape of Planet PC, it seems like a good time to recall some memories of the last major transition in the computer world. I call this the Pre-Cambrian explosion of Personal Computers.  It took place from the mid-70's to the early 80's, when dozens of entrepreneurial companies sprang up all hoping to hit the magic combination of hardware and software and sell a zillion units.

In a way, it was exactly the inverse of what is happening now. Back then, a myriad of bizarre life forms sprang up to fill the niche exposed by the advent of the microprocessor, only to be driven into extinction by the arrival of a behemoth with superior reproductive capability.

Whereas today, the new ecological zone was initially colonized by a single fierce predator, who may yet out-compete all other arrivals trying to crawl out onto the shore of the new territory.

On to a series of occasional posts on the Pre-Cambrian PC era...

Java's vision for the future

Just noticed this slideshow from QCon 2012 on the vision for Java's evolution over the next few versions.

Some highlights of JDK 8:

• closures via lambda expressions
• filter-map-reduce framework
• modules

and further out:

• JNI goes away!
• support for big in-memory data (arrays) and data structure optimizations

A key driver of these enhancements is the transition to multi-core and the end of the free performance lunch.  This is a great graph that puts it in perspective - those flat lines are the lunch bill arriving.

QCon always has excellent presentations, so I'm sure there's some other thought-provoking slideware on that site as well.

Battle of the Spatial Scripting Languages!

The gloves are down in this exchange on gis.stackexchange, over LOC counts for Fiona, JEQL, and GeoScript-JS.

Fiona takes one on the chin in the first round!  But Sean recovers quickly, and makes a hit by pointing out that Fiona allows defining new functions right in the scripting language.  This would be feasible in JEQL as well, by providing hooks to into JVM-based scripting languages.  And it's also very easy to develop extensions as simple Java functions or classes.

Now, this is really just a schoolyard scuffle over a very simple use case for spatial scripting.  What the fans really want to see is a face-off over a more complex task...

A cool application of JTS Voronoi diagrams

I just noticed this series of posts by Stephen Mather on generating medial axes of polygons by creating Voronoi diagrams on the densified linework. 

If I read correctly, he's doing this via Python GeoScript, which uses JTS for the heavy lifting.

As he says, this is a venerable technique ( hack) for creating medial axis approximations.  Regardless of its computational soundness (and hey, whatever works), it's a great example of using a series of complex geometric operations to compute a useful result.

I'll have to try this in JEQL....

(Stephen seems to have moved on to implementing the Scale Axis Transform in PGSQL - very impressive!)

Barnes Analysis for Surface interpolation

Recently I've been working on generating surfaces from irregular sets of data observations, for thematic rendering purposes.  The data I'm working with is meteorological measurements (e..g max and min temperature, average rainfall, etc.).  Here's a sample:

Maximum Daily Temperatures for November 30, 2010

There are many, many different ways of interpolating surfaces - Wikipedia has a long list.  It turns out that a classic approach for interpolation of meteorological data is Barnes Interpolation (also called Barnes Analysis).  Barnes Interpolation is a surface estimating technique which uses an initial pass to provide a smoothed estimate of the surface, and optionally further refinement passes to improve the surface to more closely match the observation values.  The initial pass produces an averaged (smoothed) surface, using a summation of exponential (Gaussian) decay functions around each observation point.  Subsequent refinement passes compute an error surface using the delta between the previous estimated surface and the observations.  The error surface is added to the previous estimate surface to refine the estimate (by reducing the delta to the observations).

I'm more of a vector guy than a raster guy, but I figured that Google and Wikipedia would make short work of this - after all, how hard can raster be?  Unfortunately, the Wikipedia article is not very helpful - the description is unclear, and it doesn't even define all the symbols used in the presented equations.  I found some better references here and here, but they still contained some confusing discrepancies in the equations.  Barnes' original paper is online as well, but as is often the case with primary sources, it is more concerned with the physical properties of the analysis, and doesn't have a very crisp description of the actual algorithm.

Triangulating between the various sources, I eventually arrived at the following algorithm definition (provided herewith without warranty!):

Barnes Interpolation Algorithm For the first pass, the estimated value Eg at each grid point gx,y is:    Eg = ∑( wi * oi) / ∑( wi ) where: oi is the value of the i'th observation point wi is the weight (decay function) value for the i'th observation point, defined as       wi = exp(- di2   /   L2C  ) where: di is the distance from the grid point to the i'th observation point L is the length scale, which is determined by the observation spacing and the natural scale of the phenomena being measured.   C is the convergence factor, which controls how much refinement takes place during each refinement step.  In the first pass the convergence factor is 1. For subsequent passes a value in the range 0.2 - 0.3 is effective. During refinement passes the estimate at each grid point is re-computed using:   E'g = Eg + ∑( wi * (oi - Ei) ) / ∑( wi ) where: Ei is the estimated value at the grid cell containing the i'th observation point

The target development platform is Java (of course).  For prototyping purposes I used JEQL to read sample datasets, run the algorithm with various parameter choices, and produce plots to visualize the results.  Here's the surface plot for the dataset above:

This looks pretty reasonable, and provides a nice visualization of the large-scale data trend. (Surprise!  It's hotter down south...)  The ultimate goal is much more exciting.  The algorithm will be provided as a Rendering Transformation in GeoServer, allowing surfaces to be generated dynamically from any supported datastore.  This will open the way to providing all kinds of further surface visualizations, including other interpolation techniques - as well as everyone's favourite, Heat Maps (aka Multivariate Kernel Density Estimation).

Now, the pure mathematical description of the algorithm is sufficient for a proof-of-concept implementation producing single images.  But to make it usable in a dynamic way in a high-demand environment, there are more details which need to be worked out to improve performance and functionality.  For instance:

• How large an extent of data points is necessary to produce a stable surface for a given area?  This is important for dynamic presentation, since zooming, panning and tiling should not change the appearance of the surface.  Of course, all data points could always be used, but more selectivity improves performance.
• How far should the surface be extrapolated into areas of few or no data observations?  Apart from the validity implications, limiting the extent of the computed surface provides a boost in performance. 
• A further performance increase may be gained by computing the surface on a coarse grid, then up-sampling in an appropriate way to the display resolution.
• For data which is global in scope, how does the algorithm need to be changed to accommodate computation on the spheroid?  It seems likely that it's necessary to use a geodetic distance calculation, and there may be other impacts as well.

Let's do it again at FOSS4G-NA!

FOSS4G 2011 in Denver last September was such a good time, we're going to do it all over again at FOSS4G North America 2012! 

It will be interesting to see who and how many attend.  Some conference budgets might have been blown last year - but note that this event is cleverly timed to fall into a new fiscal!  And being in D.C. there should be lots of suits and spooks in attendance.

I'm definitely looking forward to building on the momentum from the last conference, especially given my new gig with OpenGeo (ok, I'm definitely not booth-babe material - but being a booth-geek is probably more fun). 

Your Favourite Map Projection

...and what it says about you.

XKCD goes geospatial!  (Happens to everything eventually...)

OpenGeo Interview

My name in lights once again.  If this keeps up I'm going to have to start asking for an outrageous rider...

Sobering thought

According to this interesting BBC page on population, the number of humans on Earth has more than doubled in my lifetime.

Anyone who has spent time trying to speed up an O(n^2) algorithm knows this can't be a good thing...

Using a web service for Elevation from JEQL

MapQuest has a whole range of useful and freely available web services for geospatial queries.  One of these is the Open Elevation Service.  I thought it would be interesting to try it out using JEQL's capabilities of querying URLs and plotting datasets as images.

The Open Elevation Service is extremely simple to use.  You just submit a URL containing the location:

http://open.mapquestapi.com/elevation/v1/getElevationProfile?callback=foo&shapeFormat=raw&latLngCollection=39.740112,-104.984856

and you get back a JSON (or XML) response with the elevation:

 foo({
"shapePoints":[39.740112,-104.984856],
"elevationProfile":[{"distance":0,"height":1616}],
"info":{"copyright":{"text":"© 2011 MapQuest, Inc.","imageUrl":"http://tile21.mqcdn.com/res/mqlogo.gif","imageAltText":"© 2011 MapQuest, Inc."},
"statuscode":0,
"messages":[]}
});

JEQL doesn't (yet) do JSON parsing, but it's easy to pick out the height attribute with a simple RegEx:

height = RegEx.extract(elevJson, \'"height":(-?\d+)\}');

Querying over a grid of locations provides a dataset which is a gridded Digital Elevation Model (DEM).  JEQL provides a function to generate a table which is a grid of integers, so it's easy to produce a grid of lat/long locations:

loc = select baseLat + i * cellSize lat, baseLon + j * cellSize lon
    from Generate.grid(1, gridSize, 1, gridSize);

Once the elevations for all locations are retrieved,  the locations are used to compute square cells (polygons) which cover the area, and the cells are symbolized with a fill value from a suitable color map based on elevation.

Here's the final result.  It's slightly squashed in aspect ratio, since the data grid is in the lat/long coordinate system.  The obvious data dropouts seem to be hard errors in the MapQuest service, which is slightly surprising.  But's it's recognizably Mt. St. Helens - you can even see the lava dome in the crater!

Not surprisingly, this isn't quick to generate - the service runs at about 10 queries per second.  (Actually, that seems pretty reasonable for a free web service!).  The service does support querying multiple points per request - I might try to see if this decreases the latency (although this might be a bit complex to express in JEQL).

Here's the whole script:

//--------------------------------------
//
// Plots a DEM generated using the MapQuest Elevation service
//
//--------------------------------------
baseLat = 46.16 ;  baseLon = -122.25;    // Mt St Helens
areaSize = 0.10;
gridSize = 100;
cellSize = areaSize / gridSize;

//---- Generate location grid
loc = select     baseLat + i * cellSize lat,
        baseLon + j * cellSize lon
    from Generate.grid(1, gridSize, 1, gridSize);

//---- Query elevation service for each grid point
telev = select lat, lon, elevJson, elev
    with {
        url = $"http://open.mapquestapi.com/elevation/v1/getElevationProfile?callback=foo&shapeFormat=raw&latLngCollection=${lat},${lon}";
        elevJson = Net.readURLnoEOL(url);
        height = RegEx.extract(elevJson, \'"height":(-?\d+)\}');
        elev = Val.toInt(height);
    }
    from loc;
Mem telev;

//---- Create raster cell boxes, symbolize with color map based on elevation, and plot
minElev = 800.0; maxElev = 2500.0;
tplot = select Geom.createBoxExtent(lon, lat, cellSize, cellSize) cell,
        Color.interpolate("00a000", "ffb000", "aaaaaa", "ffffff", index) style_fillColor
    with {
        index = (elev - minElev) / (maxElev - minElev);
    }
    from telev;

textent = select geomExtent(cell) from tplot;

Plot     extent: val(textent)
    data: tplot
    file: "dem.png";

(r-i-p (john mccarthy))

John McCarthy, inventor of LISP, has died at age 84. 

Having finally balanced his parentheses, his atoms can now be returned to the free list.

Jobs, Ritchie, McCarthy - so many of the greats passing....

Extracting GTFS data using JEQL

Once you understand the General Transit Feed Specification (GTFS) data model, the next step is to transform the data into a format usable in GIS tools.  As usual, that means shapefiles (sigh...).  Here's a JEQL script that extracts the geometry of all unique trips (aka "route patterns"), as well as the stop locations.

Some interesting aspects of the script:

• the GTFS spatial data is in geographic coordinates (lat/long). The script projects this into BC-Albers (EPGS code 3005) using the CRS.project() function.  Many other coordinate systems could be used (or none - e.g. for output to KML)
• the geomConnect()aggregate function connects Points into LineStrings

/*=======================================================
Extracts Trips (with shapes) and stops from from a GTFS dataset. 
Includes shape, trip and route information.
Spatial data is converted from Geographic to BC-ALbers

Usage Notes:
- indir: directory containing the input GTFS files
========================================================*/

indir = "";  outdir = "";

//========================================================
//  Extract Trips with unique linework, with route information
//========================================================
CSVReader tshape hasColNames: file: indir+"shapes.txt";
CSVReader ttrip hasColNames: file: indir+"trips.txt";
CSVReader troutes hasColNames: file: indir+"routes.txt";  
Mem troutes;

crs_BCAlbers = "EPSG:3005";

//--- extract unique shapes
trips = select distinct route_id, trip_headsign, shape_id 
    from ttrip;

//--- extract trip pts and order them in sequence
tpt = select shape_id, Val.toInt(shape_pt_sequence) seq,
		Geom.createPoint(Val.toDouble(shape_pt_lon), 
                                 Val.toDouble(shape_pt_lat) ) pt
	from tshape
	order by shape_id, seq;

//--- connect pts into lines and project
tline = select shape_id, 
     CRS.project(Geom.toMulti(geomConnect(pt)), crs_BCAlbers) geom
	from tpt
	group by shape_id;
	
//--- join route info to trip shape
trouteLines = select r.route_id, r.agency_id, 
    r.route_short_name, r.route_long_name,
    r.route_type, r.route_url,
    trip_headsign, l.*
        from tline l
        join trips t on t.shape_id == l.shape_id
        join troutes r on r.route_id == t.route_id;

ShapefileWriter trouteLines file: outdir+"gtfs_trips.shp";

//========================================================
//  Extract Stop points
//========================================================
CSVReader tstopsRaw hasColNames: file: indir+"stops.txt";

tstop = select stop_id,stop_code,stop_name,stop_desc,
        CRS.project(Geom.createPoint(
	    Val.toDouble(stop_lon), 
	    Val.toDouble(stop_lat) ), crs_BCAlbers) pt,
        zone_id
        from tstopsRaw;

ShapefileWriter tstop file: outdir+"gtfs_stops.shp";

Using this script on the GTFS data for Metro Vancouver (available from TransLink here) produces a dataset that looks like this:

 

JTS & GEOS Wordles

For my JTS talk at FOSS4G I created a couple of Wordles showing systems which use JTS and GEOS.  Here they are for posterity (and because I love seeing Wordles...)