Contouring Data

It’s been a while since I did any Fortran. I’ve been looking into contouring algorithms and decided to have a look at Paul Bourke’s Conrec program that was originally published in Byte magazine in 1987:

http://paulbourke.net/papers/conrec

Simple Contours

The graph above shows the underlying data values as a coloured square grid with the black contour lines on top. The data point is in the centre of the grid square. Blue indicates a data value of 0.0 while red is 1.0. Contour lines are drawn for the 0.4, 0.6 and 0.8 intervals.

It is a very simple and compact algorithm, so I ended up with another C# implementation relatively quickly. There is already a C# port, along with Java, C and C++, so this was really just an aid to understanding.

Complex Contours

Contouring algorithms can be classified into one of two types: regular grids or irregular grids. The Conrec algorithm is a regular grid contour algorithm as the data values are a 2D matrix. The x and y axes can be logarithmic or irregular, but there are data values for every point on the grid.

In contrast, irregular contouring algorithms take a list of points as input and contour from them directly. This is the situation we are in with most of our GENeSIS data, but the first step in irregular grid contouring is to understand the regular grid case. The next step is to take the point data, create a Delaunay triangulation and apply the same ideas from the regular grid case, but to the triangulation.

Having looked at regular grid contouring, the next step is an implementation of Delaunay triangulation, followed by Voronoi, which is the dual of Delaunay and can be used for adjacency calculations on polygonal areas.

GMapCreator and Google API v3

I’ve created a new html template for the GMapCreator that uses Google’s API v3. In addition to not requiring the API key which is locked to the URL, this also means that you can take advantage of Google’s new styled maps base layers and Fusion table overlays.

The following map shows UK geology as the green overlay with a styled base layer called ‘Moody’:

The green data overlay shows UK geology from GMapCreator tiles while the base layer is using the 'Moody' style

Thanks to Steven Gray for the Moody style: http://bigdatatoolkit.org/

The GMapCreator template for creating Google API v3 maps can be downloaded from the following link: http://www.maptube.org/downloads/html-templates/template-APIv3.html 

(you need to right click and ‘save as’, otherwise it will open in the browser).

Select this as the html template in the GMapCreator and it will create Google API v3 maps from shapefiles automatically. I’ll post some alternative base map styles once I’ve had a chance to experiment with it some more.

Using the GMapCreator with 64 Bit Java

Using the GMapCreator on a 64 bit laptop recently, I found myself without access to a 32 bit Java Virtual Machine. As the Windows native version of the Java Advanced Imaging (JAI) project that the GMapCreator depends on is 32 bit only, I needed to use the pure Java version of JAI. It’s not immediately obvious how to do this, so I’ve detailed the method as follows:

1. Download and unpack the pure Java version of JAI 1.1.3. As the JAI project web page has recently changed, I’ve put a link to this on the MapTube website. You can download the JAI package from the following link: http://www.maptube.org/downloads/JAI/jai-1_1_3-lib.zip

2. Extract the files ‘jai_codec.jar’ and ‘jai_core.jar’ from the ‘lib’ directory and place them in the same folder as the gmapcreator.jar file e.g. C:\Program Files\CASA-UCL\GMapCreator if you are using the default installation folder on Windows.

3. Create a file called ‘run-gmc.bat’ in the same directory as GMapCreator.jar containing the following line:

[csharp]java -Xms1024M -Xmx1024M -classpath jai_core.jar;jai_codec.jar;gmapcreator.jar gmapcreator.GMapCreatorApp[/csharp]

4. Run the GMapCreator by double clicking on the ‘run-gmc.bat’ file.

The same technique can be adapted to work on Linux or Mac. The main advantage of using a 64 bit JVM is that you now have access to a 64 bit address space, so the GMapCreator isn’t limited by the 3GB limit in 32 bit any more.

Weather Underground

I’ve been looking at the Weather Underground API (http://wiki.wunderground.com/index.php/API_-_XML) which gives access to the observation stations and the data they are collecting.

All the stations returned from the Weather Underground XML API when using "London" as the search string. Colour indicates air temperature with blue=12.7C, green=13.9C and red=20.5C

The API uses simple commands to query for a list of stations, for example:

http://api.wunderground.com/auto/wui/geo/GeoLookupXML/index.xml?query=london,united+kingdom

Using C# and .net, this is accomplished as follows:
[csharp] WebRequest request = WebRequest.Create(string.Format(GeoLookupXML, @"london,united+kingdom"));
HttpWebResponse response = (HttpWebResponse)request.GetResponse();
XmlDocument doc = new XmlDocument();
doc.Load(response.GetResponseStream());[/csharp]
Then the returned XML document is parsed using XQuery to extract the station name, lat/lon coordinates and whether it is an ICAO station or a personal weather station.
[csharp]XmlNodeList Stations = doc.GetElementsByTagName("station");
foreach (XmlNode Station in Stations)
{
XmlNode IdNode = Station.SelectSingleNode("id");
XmlNode ICAONode = Station.SelectSingleNode("icao");
}[/csharp]
This gets us a list of stations ids and ICAOs which can then be used to build individual queries to obtain real time data from every station:
[csharp]foreach (string Id in PWSStations)
{
XmlDocument ob = GetCurrentPWSOb(Id);
XmlNode Ntime = ob.SelectSingleNode(@"current_observation/observation_time_rfc822");
XmlNode Nlat = ob.SelectSingleNode(@"current_observation/location/latitude");
XmlNode Nlon = ob.SelectSingleNode(@"current_observation/location/longitude");
XmlNode NairtempC = ob.SelectSingleNode(@"current_observation/temp_c");
string time = Ntime.FirstChild.Value;
string airtempC = NairtempC.FirstChild.Value;
string lat = Nlat.FirstChild.Value;
string lon = Nlon.FirstChild.Value;

//do something with the data…
}

//NOTE: only slight difference in xml format between PWS and ICAO
foreach (string ICAO in ICAOStations)
{
XmlDocument ob = GetCurrentICAO(ICAO);
XmlNode Ntime = ob.SelectSingleNode(@"current_observation/observation_time_rfc822");
XmlNode Nlat = ob.SelectSingleNode(@"current_observation/observation_location/latitude");
XmlNode Nlon = ob.SelectSingleNode(@"current_observation/observation_location/longitude");
XmlNode NairtempC = ob.SelectSingleNode(@"current_observation/temp_c");
string time = Ntime.FirstChild.Value;
string airtempC = NairtempC.FirstChild.Value;
string lat = Nlat.FirstChild.Value;
string lon = Nlon.FirstChild.Value;

//do something with the data…

}[/csharp]
After that it’s simply a matter of writing all the data to a CSV file so that you can do something with it.

Air temperature for London plotted using the MapTubeD heatmap tile renderer