Often when dealing with spatial data recorded using Ordnance Survey (OS) National Grid References (NGRs for short), you will find that different objects within your data are recorded to differing levels of spatial precision. This can lead to misrepresentative results if treated simply as point objects (especially for objects poessessing only a very vague location) and, therefore, we need a method of representation that is truer to the actual spatial precision of the data objects.
To help explain, NGRs are formatted as follows (fee free to skip this paragraph if you understand NGRs): the OS divides the UK up into 100km by 100km grid squares, which are designated using two letters. Therefore, an NGR consists of two letters followed by a string of numbers, the numbers representing the object’s location within the large grid square. There should be an even number of these digits, with the first half representing the easting and the second half the northing. However, the exact number of digits will vary according to how precisely located the object is. For example, an NGR such as SQ2222233333 is in the square SQ at the location 22222 metres east, 33333 metres north; however we would have to add zeroes to a shorter NGR as the precision of the point is less well known, e.g. SP444555 is in the square SP at the location 44400, 55500. For more on converting NGRs to numeric coordinates, see my previous post.
However, simply counting the digits is not always a sufficient measure of coordinate precision, as the person who entered the data may have tagged zeroes onto the end of the easting and northing. For example, the NGR ST5500066000 appears to be accurate to the nearest metre, but it is unlikely that a point would fall exactly onto the 55000, 66000 point and, as such, it is much more likely that this point is recorded to the nearest kilometre. I have a small Python script to perform this test, which can be downloaded here.
In any event, the end result of this is that an NGR that is apparently a point is, in fact, properly the origin point (southwest corner) of a square of width and height equal to its precision, e.g. the NGR ST5566 represents a square 1km by 1km spreading northeast from its location, whereas the NGR SQ2222233333 is a 1m by 1m square. Therefore, if we want to represent these points properly in a GIS environment, we have to present them as squares of appropriate size.
So how can we automate creating this representation in ArcMap (version 10)?* I will assume reasonable familiarity with the software.
If we have a point shapefile of our locations (use a duplicate to preserve the original), with their associated attribute data including their spatial precision**, we can create these squares in four relatively simple steps. Here are our points on a map:
First, we open up the attribute table for the layer. Then we add a new field to the table, which we shall call MIDPOINT and which will be of the Float type (i.e. a floating point number). Then we open up the Field Calculator for the MIDPOINT field and enter the calculation [Precision field] / 2.0 as follows:
This new field will be used to add on to the x and y coordinates of each point in the layer to move them to the centre of each square to be created. We then open the Field Calculator for the Shape field and enter the code snippet seen here (make sure you have the Python parser selected at the top):
The expression can be downloaded here. This results in the points moving to a new location to the north east in the midpoint of the squares we will soon create (notice that some points have moved a long way yet others remain close to their origin):
Next we create buffers for our moved points using the Buffer tool in ArcToolbox. You want to set it to use the MIDPOINT field as the buffer size (we use the MIDPOINT rather than the full precision field as the size will be the radius of the circle created, not the diameter) and you do not want it to dissolve the buffers where they overlap, as we want one object for each point. The result should look something like this:
Finally, we use the Feature Envelope to Polygon tool in ArcToolbox using the buffer layer as input. This, helpfully, retrieves the envelopes that enclose each circle (used to make spatial queries faster etc.), which is equivalent to the squares that we were after in the first place. Therefore, this layer is the result that we desired:
Any points that do not seem to have output squares in the image do, in fact, have them but they are too small to see at this scale (e.g. 1m or 10m precision). Overall, this is a simple way to do this task in just a few minutes whilst avoiding any overly complex scripting. I think I will work on a Python script to achieve this task in one step, but the method just described works perfectly well in the interim.
The final result means that we now both have a visual representation of the area to which each of our NGR “points” actually relates, and also a layer which can be used in further analyses when trying to understand the contents of our layer: all of the layer attributes from the original layer should find their way automatically into the final output, so we can perform whatever analyses we wish on our new, more properly representative geographic objects.
* ArcMap is part of ESRI’s ArcGIS and is the GIS software that we are predominantly using in this project.
** If you did not record the precision of the NGRs whilst converting them, you could use the Field Calculator to work it out (in the same way as described above for the other calculations). You would have to create a fields called something like ‘PRECISION’ and ‘TRUEPRECIS’, and then use the expression code downloadable here and here in turn on the first and then second fields.