Category: GIS

Some time ago, I looked at some uses for Ordnance Survey Open Data, coming to the conclusion that a sensible way to work with it would be to import it in to a geospatial-enabled database.
Each set of data is provided in ESRI Shapefile format, and has four files:

• shp – shape format
• shx – shape index format
• dbf – attribute format in dBase IV format
• prj – projection format

The shp2pgsql command converts SHP files in to a set of SQL commands which will effectively import the data in to PostgreSQL. Here’s a ridiculously simple guide to importing a file:
createdb os_opendata
psql -d os_opendata -c "CREATE EXTENSION POSTGIS"
shp2pgsql <filename>.shp <table_name> | psql -d os_opendata
Depending on the speed of your machine, in a few seconds you’ll find a new table in your database with all the data included.
And finally, what if you just want to import all of the data at once? Try this:
find Data/ -name "*.shp" | xargs -I % -n1 shp2pgsql % | psql -d os_opendata

The excellent people at the Greater London Assembly have released a list of bus stops and bus routes in London. The coordinates of each bus stop are in eastings and northings, and I wanted to convert these to longitude and latitude for my Ruby on Rails application.
Using the proj4rb gem and some projection definitions from spatialreference.org – with some help from Harry Wood’s blog, I came up with the following code:

#!/usr/bin/ruby
require 'rubygems'
gem 'proj4rb'
require 'proj4'
easting = 529978
northing = 186491
srcPoint = Proj4::Point.new(easting, northing)
srcProj = Proj4::Projection.new('+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs')
dstProj = Proj4::Projection.new('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
dstPoint = srcProj.transform(dstProj, srcPoint)
print "lat=#{dstPoint.lat * (180 / Math::PI)}n"
print "lon=#{dstPoint.lon * (180 / Math::PI)}n"

To convert WGS84 to OSGB36:

#!/usr/bin/ruby
require 'rubygems'
gem 'proj4rb'
require 'proj4'
latitude = 51.5623279577278
longitude = -0.126277004538848
srcPoint = Proj4::Point.new(longitude * (Math::PI / 180), latitude * (Math::PI / 180))
srcProj = Proj4::Projection.new('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
dstProj = Proj4::Projection.new('+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs')
dstPoint = srcProj.transform(dstProj, srcPoint)
print "easting=#{dstPoint.x}n"
print "northing=#{dstPoint.y}n"