Archives for posts with tag: Python

I recently discovered the GeoAlchemy2 project – a replacement for the original GeoAlchemy package, focused on providing PostGIS support for SQLAlchemy. The SQLAlchemy package is a “Python SQL Toolkit and Object Relational Mapper”. In a nutshell this means you can write Python classes and map them to PostgreSQL tables without the need to write SQL statements – pretty cool!

PostGIS is great for doing spatial stuff, but if you’re using it as back-end for a Python app then you can spend a lot of time writing Python wrappers around SQL statements, and even with the excellent Psycopg2 package this can be tricky. This is especially true if you’re using the OGR Python bindings to handle PostGIS read/writes.

Enter GeoAlchmey2. I’ve been experimenting with it for a week, in that time I’ve learnt this:

For developing geospatial Python apps with PostGIS, GeoAlchemy2 is nothing short of revolutionary.

You can call PostGIS functions in Python, which means you can use them (and the data) directly within your Python application logic. Here’s an example. The SQL statement below uses PostGIS to create a new line geometry between a point, and the closest point on the nearest line.

Now here’s a snippet from a Python script, performing the same process using GeoAlchemy2.

You’ll notice here that we’re actually calling our own Python function “make_link_line” during the query to create the new geometry. This exemplifies how we can move PostGIS objects around inside the script. Once the query runs we can access the returned data in our application from the row variable. Below is the complete script.

Nearest neighbour PostGIS and GeoAlchemy2 script:

The script above is just a simple example, but it shows how powerful GeoAlchemy2 is for embedding PostGIS objects and methods inside Python. I’m really looking forward to digging deeper into the functionality of GeoAlchemy2 and SQLAlchemy to integrate them within my own projects. Check out the official tutorials for more examples:


Plotting confidence intervals of linear regression in Python

After a friendly tweet from @tomstafford who mentioned that this script was useful I’ve re-posted it here in preparation for the removal of my Newcastle University pages.

This script calculates and plots confidence intervals around a linear regression based on new observations. After I couldn’t find anything similar on the internet I developed my own implementation based on Statistics in Geography by David Ebdon (ISBN: 978-0631136880).

Linear regression plot

Plot of linear regression with confidence intervals

# - example of confidence limit calculation for linear regression fitting.

# References:
# - Statistics in Geography by David Ebdon (ISBN: 978-0631136880)
# - Reliability Engineering Resource Website:
# -
# - University of Glascow, Department of Statistics:
# -

import numpy as np
import matplotlib.pyplot as plt

# example data
x = np.array([4.0,2.5,3.2,5.8,7.4,4.4,8.3,8.5])
y = np.array([2.1,4.0,1.5,6.3,5.0,5.8,8.1,7.1])

# fit a curve to the data using a least squares 1st order polynomial fit
z = np.polyfit(x,y,1)
p = np.poly1d(z)
fit = p(x)

# get the coordinates for the fit curve
c_y = [np.min(fit),np.max(fit)]
c_x = [np.min(x),np.max(x)]

# predict y values of origional data using the fit
p_y = z[0] * x + z[1]

# calculate the y-error (residuals)
y_err = y -p_y

# create series of new test x-values to predict for
p_x = np.arange(np.min(x),np.max(x)+1,1)

# now calculate confidence intervals for new test x-series
mean_x = np.mean(x)			# mean of x
n = len(x)				# number of samples in origional fit
t = 2.31				# appropriate t value (where n=9, two tailed 95%)
s_err = np.sum(np.power(y_err,2))	# sum of the squares of the residuals

confs = t * np.sqrt((s_err/(n-2))*(1.0/n + (np.power((p_x-mean_x),2)/

# now predict y based on test x-values
p_y = z[0]*p_x+z[0]

# get lower and upper confidence limits based on predicted y and confidence intervals
lower = p_y - abs(confs)
upper = p_y + abs(confs)

# set-up the plot
plt.xlabel('X values')
plt.ylabel('Y values')
plt.title('Linear regression and confidence limits')

# plot sample data
plt.plot(x,y,'bo',label='Sample observations')

# plot line of best fit
plt.plot(c_x,c_y,'r-',label='Regression line')

# plot confidence limits
plt.plot(p_x,lower,'b--',label='Lower confidence limit (95%)')
plt.plot(p_x,upper,'b--',label='Upper confidence limit (95%)')

# set coordinate limits

# configure legend
leg = plt.gca().get_legend()
ltext = leg.get_texts()
plt.setp(ltext, fontsize=10)

# show the plot

Visualising London’s buildings with 3D mapping

I recently discovered a JavaScript library called OSM Buildings which provides a 3D representation of buildings from Open Street Map (OSM) data on a 2D map. OSM Buildings doesn’t generate a complete 3D model but uses a clever shading technique to make buildings appear ‘3D’. OSM Buildings and can be used with either OpenLayers or the Leaflet JavaScript map API’s.

3DLondon map screenshot

I was interested in whether OSM Buildings could be used to visualise non-OSM data, in particular the Geoinformation group’s UKMap product, to create a complete ‘3D’ city map. The UK map database is a commercial, attribute rich, map database designed to compete with the UK’s Ordnance Survey products, and contains a LiDAR derived height for most structures in urban areas. Furthermore, UK academics can access the UKMap data via the MIMAS Landmap service. As part of my PhD research I obtained a copy of the UKMap basemap and building heights for London from the Landmap team.

The OSM Buildings API provides a method to add a buildings layer from a GeoJSON object, so I wrote a simple Python CGI script to grab building features from the UKMap PostGIS database and return it as GeoJSON. The embedded SQL query (via the excellent Psycopg module) uses the new PostGIS 2.0 function ST_AsGeoJSON to return a JSON representation of the geometry. Unfortunately this isn’t as straight forward as it seems as PostGIS doesn’t return any feature attributes. However after reading this blog at the Postgres online journal and installing Andrew Dunstan’s backport of JSON support for PostgreSQL 9.1, I managed to get the query to return proper JSON representation of the attributes and geometry.

To reduce data volume (there are ~5 million features in the database) the CGI script accepts a bounding box parameter in well known text form, and only features within the bounding box are returned. Here’s a snippet showing the SQL query:

SELECT row_to_json(fc)
FROM (SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
      FROM (SELECT 'Feature' As type,
            ST_AsGeoJSON(ST_Transform(lg.the_geom,4326))::json As geometry,
            row_to_json((SELECT l FROM (SELECT height) As l)) As properties
            FROM "3DLondon" As lg
            WHERE ST_Within(lg.the_geom,
            LIMIT 100000)
      As f )
As fc;',[bbox_wkt])

A sample of the features and their geometry returned as GeoJSON:


Thanks to some good documentation at OSM Buildings the front end is a simple JavaScript function to add the buildings to the map via a jQuery post AJAX call.

//Request data from GeoJson server and add to map
function getBuildings(bounds_wkt){
	        .setStyle({ wallColor: 'rgb(200,190,180)', roofColor: 'null',
                        strokeRoofs: true })

Using this function it was fairly straightforward to add some functions to listen for Leaflet map events to refresh the building data when the map was moved or zoomed.

//Listen for map pan events and get new buildings data
map.on('moveend', function(e) {

From the screenshot above you can see that OSM Buildings is doing a great job of rendering stuff in ‘3D’. The basemap is Open Street Map data, styled and served by CloudMade with OSM building footprints removed for clarity. In a nutshell I was pretty impressed with both the OSM Buildings and the Leaflet libraries. You can check out the complete source code at, and here’s another screenshot over a larger area with the buildings rendered in blue.

3DLondon map screenshot

3DLondon map of St Paul’s/The City and surrounds

Following in the footsteps of the Raster Processing Suite, I’ve added GitHub pages for the PostGISDroid script which I wrote last year. This was a prototype Python script, built on the Android scripting layer, to log the location of an Android device to a remote PostGIS server. Check it out here:

PostGISDroid track in QGIS
A track from PostGISDroid in Quantum GIS.

If you’re looking for a full-blown Android client for mobile data acquisition check out the Sense Cloud Framework:

As mentioned previously I’m in the process of migrating content from my Newcastle pages to this site. Given that I’m already hosting most of my software projects on GitHub, I’ve taken the leap and created GitHub pages for the Raster Processing Suite QGIS plug-in. The Raster Processing Suite is a nifty bit of software I developed as part of my PhD thesis, to help with exploratory numerical analysis of time series AVHRR satellite imagery (see Holderness et al., 2013). Part of the PyRaster tools (more on those soon) created during my PhD, the Raster Processing Suite is a raster calculator with a Python scripting interface (including automatic script generation) for raster image processing in Quantum GIS (see screen-shot below).

GitHub pages provide free web hosting for content related to software projects. You can check out the Raster Processing Suite at Furthermore, with a bit of tweaking I’ve also managed to get the pages to host the QGIS XML repository for direct download into QGIS:

I’ll post links to GitHub pages for other projects in the coming weeks.

Screenshot of the Raster Processing Suite QGIS plug-in

Screen-shot of the Raster Processing Suite QGIS plug-in