Tags: gis
Calculating distance with Mysql
By admin on Sep 30, 2009 | In Formula 1 - More advanced stuff, Showroom - Examples of Mysql usage | Send feedback »
After the recent posting of data in Wikileaks, it's claimed that here in the UK we can now freely access our longitude and latitude co-ordinates for each post code. However, some of the data in their files may well be deliberately corrupted - you may be infringing copyright law by using that data, and the mistakes will prove that you've used illegal data. Free the Postcode is perhaps a better place to get free data - http://www.freethepostcode.org/
Now, it's possible to calculate distances between each point using the following calculation (originally from the fine book Pro Mysql):
The distance d between two points (x1,y1) and (x2,y2) can be calculated from the following equation (so long as x values are latitude and y values are longitude in radians, r is the radius of the sphere which is 3956 in miles):
Code:
d= acos(sin(x1)sin(x2)+cos(x1)cos(x2)*cos(y2-y1)) * r |
Initially I downloaded relevant postcode data file onto my Linux system. I then (by using regex in egrep) selected those post codes beginning NG and a single digit (Nottingham - I'm not interested in loading every single postcode at this stage) and use the shell command cut to extract only those fields I'm after. In this case it's a comma-delimited file so we use the -d switch to let the function cut know that commas are the delimiter.
Code:
egrep 'NG[0-9]+ [0-9]' freepostcodesfile | cut -d , -f 1,14,15 > ng.txt |
I did the same for OX postcodes so I had some data in two different towns to play around with.
Prior to inserting the postcode data I used the following script to create the table (I'll explain the last three columns a little later).
Code:
CREATE TABLE postcodes ( postcode CHAR(8) NOT NULL | |
, longitude DECIMAL(9,6) NOT NULL | |
, latitude DECIMAL(9,6) NOT NULL | |
, rad_xaxis DECIMAL(9,6) | |
, rad_yaxis DECIMAL(9,6) | |
, rad_zaxis DECIMAL(9,6) | |
, PRIMARY KEY(postcode)); |
To insert the data into my table I used :
LOAD DATA LOCAL INFILE '/home/usname/gis/ng.txt' INTO TABLE postcodes FIELDS TERMINATED BY ',' OPTIONALLY ENCLOSED BY '"' LINES TERMINATED BY '\n' (postcode, longitude, latitude);
and again for the ox file (Oxford postcodes)
LOAD DATA LOCAL INFILE '/home/usname/gis/ox.txt' INTO TABLE postcodes FIELDS TERMINATED BY ',' OPTIONALLY ENCLOSED BY '"'
LINES TERMINATED BY '\n' (postcode, longitude, latitude);
After putting our data into a table we can then use the trig functions available in Mysql to calculate distances between two points. However we have to be aware of time-consuming functions, that make the selection of records from the table where the distance is below a certain figure (ie after the calculation of all relevant records) can be very cumbersome. It's better to precompute and store much of this in the database table where we're recording the postcodes. To start with we need to change the latitude and longitude figures given from degrees into radians. This is done using the following formula:
radians= degrees*(PI/180)
Luckily we can obtain PI in Mysql by using PI()
What I will do here is precompute and populate 3 extra fields for each record (in radians) and store them as rad_xaxis, rad_yaxis and rad_zaxis, as per the following definition:
xaxis = cos(radians(Lat)) * cos(radians(Lon))
yaxis = cos(radians(Lat)) * sin(radians(Lon))
zaxis = sin(radians(Lat))
So we update our calculated fields as per the following:
Code:
update postcodes set rad_xaxis | |
= cos((latitude * (PI() /180))) * cos((longitude * (PI() /180))); | |
update postcodes set rad_yaxis | |
= cos((latitude * (PI() /180))) * sin((longitude * (PI() /180))); | |
update postcodes set rad_zaxis | |
= sin((latitude * (PI() /180))); |
The distance between points can then be calculated using SQL loosely like (acos( xaxis * $xaxis + yaxis * $yaxis + zaxis * $zaxis ) * 3956) (where those starting with a $ are precomputed for the start point in question in the same manner as above)
So to get my distance between say the post code OX1 1AB (central Oxford) and those postcodes starting 'OX4 3Y%'
Code:
select b.postcode, acos(b.rad_xaxis * a.rad_xaxis + b.rad_yaxis * a.rad_yaxis + b.rad_zaxis * a.rad_zaxis) * 3956 AS distance | |
from postcodes a, postcodes b | |
where a.postcode = 'OX1 1AB' | |
and b.postcode like 'OX4 3Y%' | |
and b.latitude > 0; |
Note I don't want to calculate where the file concerned didn't have data, so add to the where clause the restriction that the latitude figure is greater than 0.
I have done some calculations using this formula (for example where the fixed post code is NG7 4ER in Nottingham), and it seems to be more accurate over greater distances - I'm not sure that you could trust it for under 5 miles. Even so it's very responsive and you can add different parts of the UK and have a distance calculator for different areas easily available. It should be remembered that as the Earth is not a perfect sphere, our calculations here will not be totally accurate, but we seem to be within 200 yards for the distance between Oxford and Nottingham (as the crow flies).
Now it's often the case that you won't want all the distances between lots of random postcodes (even within your home city). However, suppose you want to bring back all supermarkets (or takeaway shops or pubs etc) within 6 miles of your designated postcode? So long as you have the supermarkets recorded in a table along with the correct post code (make sure the format is exactly the same as the postcode table) then you can use the following statement:
Code:
select b.postcode, acos(b.rad_xaxis * a.rad_xaxis | |
+ b.rad_yaxis * a.rad_yaxis | |
+ b.rad_zaxis * a.rad_zaxis) * 3956 AS distance, | |
s.name, s.address | |
from postcodes a, postcodes b, shops s | |
where b.postcode = s.postcode and a.postcode = 'OX1 1AB' | |
and b.latitude > 0 | |
and s.shoptype = 'Supermarket' | |
and acos(b.rad_xaxis *a.rad_xaxis | |
+ b.rad_yaxis * a.rad_yaxis | |
+ b.rad_zaxis * a.rad_zaxis ) * 3956 <= 6 order by distance \G |
This brought back the following records (those records I've recorded in my shops database table - in real life obviously there are others within the specified distance):
******* 1. row *******
postcode: OX2 7BY
distance: 2.66911298653386
name: Midcounties Co-Op
address: 228-240 Banbury Road, Oxford
******* 2. row *******
postcode: OX4 3XQ
distance: 3.50779225250805
name: Midcounties Co-Op
address: 21 Templars Square, Cowley
******* 3. row *******
postcode: OX3 8RA
distance: 3.88607301210902
name: Midcounties Co-Op
address: Atkyns Road, Headington
******* 4. row *******
postcode: OX4 6XJ
distance: 4.67029242514653
name: Tesco
address: Oxford Retail Park, Cowley
4 rows in set (0.01 sec)
How good is that? I'll have to round the figures though - it doesn't make any sense to calculate to 15 decimal points! Hope you get as much enjoyment from using this technique as I do. Until next time - Enjoy your Mysql coding!