Demystifying Coordinate Transformations

Introduction

This blog is going to look at the mathematical calculations involved in rendering 2D GIS vector data, based on a particular area that we want to see on a map. I.e. converting vertices from map space into view space, and back.

For each vertex in a feature geometry, we need to perform a set of operations to calculate the screen position, based on our chosen location, scale, and rotation. Also, most applications will use top left as the origin rather than bottom left for map projections.

Normally, we would use client software such as OpenLayers, or application software such as GeoServer to make these calculations (and render the points, lines, or polygons) for us, however, occasionally it is more convenient to make the calculations required ourselves. E.g. generation of an image file on a server when the data is not available as a WMS through GeoServer.

At thinkWhere, we have implemented these calculations in a number of programming languages over the years as the languages we have used have changed. We have gone from C++, through C# and Actionscript (Flash), to Python and Javascript.

The Maths

To make the calculations we use Affine Transformations. For 2D data, these take the form of 3×3 Matrices which can be applied to a point location, producing a modified point.

We can create a number of transformations, each of which represents an action on the data. E.g translate, scale, or rotate the point.

These can be combined through matrix multiplication into a single transformation.

We can also calculate the inverse of a matrix, making it just as easy to perform the reverse set of actions.

The beauty of this approach is that we can apply a number of actions to each vertex in a small number of arithmetic operations making them quick to apply. Even for rotation, the sin/cosine values are calculated once, and the calculated matrix values used in simple multiplications and additions from then on. Once the completed matrix has been generated, the calculations required for each vertex are:

x' = (a * x) + (b * y) + c

y' = (d * x) + (e * y) + f

The Good News

The good news is that once we have created a class to perform the matrix arithmetic, we don’t have to think about it again. We can then create a sequence of actions which will transform our map coordinates to view coordinates.

  • Translate the data based on our chosen centre point – centring the data at 0,0
  • Scale the data based on our chosen view scale
  • Rotate the data around the centre if required
  • Flip the horizontal axis to switch the origin from bottom left to top left
  • Translate again to the centre of our view (based on the view size)

The example below shows this in action.

Application

Sample Application

This sample application has been written using an HTML5 Canvas to draw a GIS dataset from a GeoJSON file. It uses the turf library to assist in reading the GeoJSON files.

It uses an implementation of affine transformations to calculate the necessary transform to convert map coordinates in the GeoJSON file to screen coordinates. A standard class called CanvasRenderingContext2D draws the features. We have added additional members to make the drawing easier. This object already implements affine transformations. You can apply individual actions using the rotate(), scale(), and translate() methods, which alter the current transformation to allow you to build up a set of actions, or you can supply the six relevant parameters of the matrix in a single call to transform().

We have used the approach of calculating the matrix ourselves, in order to also produce the inverse transformation. This is then used to make the reverse calculation I.e. for a point clicked on the map, find the map coordinate that this relates to in order to then search for the clicked feature within the source data.

The code below uses the view control values to create the required transformation, and pass this to the CanvasRenderingContext2D

function setupTransform(canvas, context)
 {
 var r = document.getElementById("rangeAngle");
 var tx = document.getElementById("rangeCentreX");
 var ty = document.getElementById("rangeCentreY");
 var sc = document.getElementById("rangeScale");
 var scalefactor = canvas.width/sc.value;

var t = new Transform();
 t.apply(new Translate(-parseFloat(tx.value),-parseFloat(ty.value)));
 t.apply(new HorizontalFlip());
 t.apply(new Rotate(parseFloat(r.value)));
 t.apply(new Scale(scalefactor));
 t.apply(new Translate(250, 250));

currentTransform = t;
 currentInverse = t.inverse();
 currentScaleFactor = scalefactor;

context.setTransform(t.matrix[0][0],t.matrix[1][0],t.matrix[0][1],t.matrix[1][1],t.matrix[0][2],t.matrix[1][2]);
/*
alternative method using a combination of the built-in methods
 context.translate(250,250);
 context.scale(scalefactor,-scalefactor);
 context.rotate(parseFloat(r.value));
 context.translate(-parseFloat(tx.value),-parseFloat(ty.value));
*/
 }

When a data file is loaded, the maximum extent of the features is calculated, and the controls are set up accordingly. When the controls are modified the transformation is adjusted, and the features redrawn.

Loops and more

Since it it quite simple to set up a set of actions passing in the various parameters, we can make use of loops to run though a variety of settings and start having some fun.

This example uses the Python libraries: Shapely (to read the GeoJSON data), Pillow (drawing to image), and images2gif (animated GIF) to produce an animated map.

Address Search OS OpenNames with PostGIS, SQLAlchemy and Python – PART 1

In this two part post we’ll look at implementing an address search using the Ordnance Survey Open Names dataset. We’ll use the power of Postgres with the PostGIS extension leveraging it’s built in Full Text Search, and use Python and the SQLAlchemy ORM to create a simple CLI.

address-search

Part 1 – Data Load and DB Config

Address Data

The UK is badly served for free address data. The best we have is the Ordnance Survey OpenNames dataset. It will work as a Postcode lookup or a street finder (at a push), but the dataset would require a lot of additional processing to be a useful address search. OS really want you to purchase AddressBase.

That said, OpenNames will suffice for this example and it should be easy to extend the example to a fuller dataset if you’re lucky enough to have one.

Loading Data to PostGIS

You can download OpenNames as either CSV, or GML. I’d recommend GML as it’s simpler to load it into PostGIS using OGR2OGR.

Once you unzip the archive you’ll see that the files are referenced according to the British National Grid, so you can load as much or as little as you want.

We’ll load NS68 which contains addresses in my home town of Stirling, as follows (swap out the values for your db):

ogr2ogr -f PostgreSQL PG:"host=localhost dbname=Real-World port=5432 user=iain password=password" NS68.gml -progress -nln open_names --config PG_USE_COPY YES 

You should now have a new table called open_names containing the addressing info.

Note if you want to load more gml files just use the -append flag:

ogr2ogr -f PostgreSQL PG:"host=localhost dbname=Real-World port=5432 user=iain password=password" NS88.gml -append -progress -nln open_names --config PG_USE_COPY YES 

Setting up Full Text Search

We now have our open_names table, but no text search column. So we can add a textsearchable column which must be of type TSVECTOR as follows:

ALTER TABLE open_names ADD COLUMN textsearchable TSVECTOR;

We can populate the column by using the built in function TO_TSVECTOR, this tokenises the words based on the supplied config, in our case english. However, multiple configs are supported.

UPDATE open_names SET textsearchable = TO_TSVECTOR('english', text || ' ' || localid);

If you look at the data in your new column you’ll see that it now contains text tokens representing the address data.

Increase accuracy by concatenating multiple columns

Note that we’re concatenating 2 columns together in this update statement – text and localid. In our case the reason for doing this is that the postcode in the localid column is stored without a space, meaning our search will return a result if the user enters a postcode without a space.

However, it should be clear if we had better address data, we could concat multiple columns. Meaning if a user searched for “1 Main St, Stirling, FK3 4GG” we would be able to return an accurate match.

Add an Index for faster searching

Now that we have data set up we can add an index to our new column which will ensure searches are fast:

CREATE INDEX textsearch_idx ON open_names USING GIN (textsearchable);

Let’s do some searches

Now lets query our new column to see if we can find some matches using the TO_TSQUERY function

SELECT COUNT(1) FROM open_names WHERE textsearchable @@ TO_TSQUERY('english', 'avenue')

Here we find we have 41 streets in Stirling area containing the word avenue. You’ll note that I don’t need to worry about lowercase, uppercase or where the word might appear in the string. Full text search takes care of that for me 🙂

The @@ operator basically means that the query matches the tsvector column.

Using AND and OR for better matches

A very powerful feature of Postgres’ Full Text Search is the ability to find matches contain all or some of the words in the query using the AND & operator or the OR | operator, as these examples show:

select * from open_names where textsearchable @@ to_tsquery('english', 'forth & view');

Here we only return one result Forth View which contains both Forth and View, if we change this to an OR search:

select * from open_names where textsearchable @@ to_tsquery('english', 'forth | view')

We get 7 results including Forth View, Bruce View, Forth Place.

Again it should be easy to see how powerful text searches could be built for complex text documents.

A final note on Triggers

While our address data is fairly static, if you had a table where users were regularly editing address data, or any other columns you wanted to run a full text search on, you should consider adding a trigger to keep the TSVECTOR column up to date, as outlined here.

So for our example the trigger would look like:

CREATE TRIGGER tsvectorupdate BEFORE INSERT OR UPDATE
ON open_names FOR EACH ROW EXECUTE PROCEDURE
tsvector_update_trigger(textsearchable, 'pg_catalog.english', localid, text);

Up Next

Hopefully Part 1 has demonstrated how it is very easy to set up powerful text searching in Postgres. In Part 2 we’ll look at how we can use Python and SQLAlchemy to allow you to integrate this functionality into your apps and APIs.