Olympic data
import json
fname = '/home/h02/itpe/support/olympic_path.json'
with open(fname) as fh:
r = json.load(fh)
# Print the first 3 visits, so that we have an idea of how the data is stored.
r['visits'][:3]
[{u'Latitude': 50.06675, u'Longitude': -5.71156, u'day': 1, u'entityId': 1000001, u'placeName': u"Land's End", u'sequence': 1, u'type': 1, u'url': u'lands-end'}, {u'Latitude': 50.07067, u'Longitude': -5.69666, u'day': 1, u'entityId': 1000002, u'placeName': u'Sennen', u'sequence': 2, u'type': 3, u'url': u''}, {u'Latitude': 50.10626, u'Longitude': -5.55002, u'day': 1, u'entityId': 1000003, u'placeName': u'Newlyn', u'sequence': 3, u'type': 3, u'url': u''}]
Convert the data into shapely geometries. Let's create a list of shapely Point objects as well as a dictionary containing a list of point objects for each day.
import shapely.geometry as sgeom
from collections import defaultdict
points = []
points_by_day = defaultdict(list)
for visit in r['visits']:
# Create the actual point.
point = sgeom.Point([visit['Longitude'], visit['Latitude']])
# Add the point to the list and dictionary.
points.append(point)
points_by_day[visit['day']].append(point)
# Make a buffered polygon of the points - this is the easiest way to visualise the data in one go.
points_buffered = sgeom.MultiPoint(points).buffer(0.05)
# Calculate the extents of the data so that we can pass that to the underlying map.
x0, y0, x1, y1 = points_buffered.bounds
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines('10m')
ax.set_extent([x0, x1, y0, y1])
ax.add_geometries(points_buffered, ccrs.PlateCarree(), edgecolor='black')
plt.show()
Wikipedia has a similar plot, suggesting that we've loaded the data correctly:
Let's create a LineString for each day's locations, and put them all into a MultiLineString.
day_routes = []
for points in points_by_day.values():
if len(points) > 1:
ls = sgeom.LineString(sgeom.MultiPoint(points))
else:
# We only had one location, so just draw a buffer around that
# one point, and get its exterior LineString.
ls = points[0].buffer(0.001).exterior
day_routes.append(ls)
all_routes = sgeom.MultiLineString(day_routes)
Let's look at these routes on a map (buffering by 0.2 degrees so that we can see it clearly):
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Mercator())
ax.coastlines('10m')
ax.set_extent([x0, x1, y0, y1])
ax.add_geometries(all_routes.buffer(0.1), ccrs.PlateCarree(),
facecolor='blue', edgecolor='black', alpha=0.2)
plt.show()
Let's do something with this data - using the Natural Earth datasets we can get hold of the UK counties:
import cartopy.io.shapereader as shpreader
counties_shp = shpreader.natural_earth(resolution='10m',
category='cultural',
name='admin_1_states_provinces')
uk_counties = []
for county in shpreader.Reader(counties_shp).records():
if county.attributes['admin'] == 'United Kingdom':
uk_counties.append(county)
# Print all UK counties sorted alphabetically.
print ', '.join(sorted(county.attributes['name'] for county in uk_counties))
Aberdeen, Aberdeenshire, Anglesey, Angus, Antrim, Ards, Argyll and Bute, Armagh, Ballymena, Ballymoney, Banbridge, Barking and Dagenham, Barnet, Barnsley, Bath and North East Somerset, Bedford, Belfast, Bexley, Birmingham, Blackburn with Darwen, Blackpool, Blaenau Gwent, Bolton, Bournemouth, Bracknell Forest, Bradford, Brent, Bridgend, Brighton and Hove, Bristol, Bromley, Buckinghamshire, Bury, Caerphilly, Calderdale, Cambridgeshire, Camden, Cardiff, Carmarthenshire, Carrickfergus, Castlereagh, Central Bedfordshire, Ceredigion, Cheshire East, Cheshire West and Chester, City, Clackmannanshire, Coleraine, Conwy, Cookstown, Cornwall, Coventry, Craigavon, Croydon, Cumbria, Darlington, Denbighshire, Derby, Derbyshire, Derry, Devon, Doncaster, Dorset, Down, Dudley, Dumfries and Galloway, Dundee, Dungannon, Durham, Ealing, East Ayrshire, East Dunbartonshire, East Lothian, East Renfrewshire, East Riding of Yorkshire, East Sussex, Edinburgh, Eilean Siar, Enfield, Essex, Falkirk, Fermanagh, Fife, Flintshire, Gateshead, Glasgow, Gloucestershire, Greenwich, Gwynedd, Hackney, Halton, Halton, Hammersmith and Fulham, Hampshire, Haringey, Harrow, Hartlepool, Havering, Herefordshire, Hertfordshire, Highland, Hillingdon, Hounslow, Inverclyde, Isle of Wight, Isles of Scilly, Islington, Kensington and Chelsea, Kent, Kingston upon Hull, Kingston upon Thames, Kirklees, Knowsley, Lambeth, Lancashire, Larne, Leeds, Leicester, Leicestershire, Lewisham, Limavady, Lincolnshire, Lisburn, Liverpool, Luton, Magherafelt, Manchester, Medway, Merseyside, Merthyr Tydfil, Merton, Middlesbrough, Midlothian, Milton Keynes, Monmouthshire, Moray, Moyle, Neath Port Talbot, Newcastle upon Tyne, Newham, Newport, Newry and Mourne, Newtownabbey, Norfolk, North Ayshire, North Down, North East Lincolnshire, North Lanarkshire, North Lincolnshire, North Somerset, North Tyneside, North Yorkshire, Northamptonshire, Northumberland, Nottingham, Nottinghamshire, Oldham, Omagh, Orkney, Oxfordshire, Pembrokeshire, Perthshire and Kinross, Peterborough, Plymouth, Poole, Portsmouth, Powys, Reading, Redbridge, Redcar and Cleveland, Renfrewshire, Rhondda, Cynon, Taff, Richmond upon Thames, Rochdale, Rotherham, Royal Borough of Windsor and Maidenhead, Rutland, Salford, Sandwell, Scottish Borders, Sefton, Sheffield, Shetland Islands, Shropshire, Slough, Solihull, Somerset, South Ayrshire, South Gloucestershire, South Lanarkshire, South Tyneside, Southampton, Southend-on-Sea, Southwark, Staffordshire, Stirling, Stockport, Stockton-on-Tees, Stoke-on-Trent, Strabane, Suffolk, Sunderland, Surrey, Sutton, Swansea, Swindon, Tameside, Telford and Wrekin, Thurrock, Torbay, Torfaen, Tower Hamlets, Trafford, Vale of Glamorgan, Wakefield, Walsall, Waltham Forest, Wandsworth, Warrington, Warwickshire, West Berkshire, West Dunbartonshire, West Lothian, West Sussex, Westminster, Wigan, Wiltshire, Wokingham, Wolverhampton, Worcestershire, Wrexham, York
Let's plot these counties, to double check that we've got them correctly located:
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.OSGB())
ax.coastlines('10m')
ax.add_geometries([county.geometry for county in uk_counties],
ccrs.PlateCarree(), facecolor='none', edgecolor='black')
plt.show()
Ok, now lets find out which counties were unlucky and did not get a visit from the Olympic torch (according to our line string):
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.OSGB())
ax.coastlines('10m')
for county in uk_counties:
if not county.geometry.intersects(all_routes):
print county.attributes['name']
# Draw the county buffered and with a low alpha first to draw the attention to the area.
ax.add_geometries([county.geometry.buffer(0.15)],
ccrs.PlateCarree(), facecolor='red', edgecolor='black', alpha=0.1)
ax.add_geometries([county.geometry],
ccrs.PlateCarree(), facecolor='coral', edgecolor='none')
ax.add_geometries(all_routes,
ccrs.PlateCarree(), facecolor='none', edgecolor='blue', alpha=0.5)
plt.show()
East Renfrewshire South Lanarkshire Trafford Isles of Scilly
And now, lets find out if there were any lucky counties which had the torch all day:
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.OSGB())
ax.coastlines('10m')
for day_route in day_routes:
for county in uk_counties:
if county.geometry.contains(sgeom.MultiPoint(day_route.coords[:])):
print county.attributes['name']
# Add an outline of the geometry with a low alpha value to emphasise
# the area.
ax.add_geometries([county.geometry.buffer(0.15)],
ccrs.PlateCarree(), facecolor='red', edgecolor='black', alpha=0.1)
ax.add_geometries([county.geometry],
ccrs.PlateCarree(), facecolor='coral', edgecolor='none')
ax.add_geometries(all_routes,
ccrs.PlateCarree(), facecolor='none', edgecolor='blue', alpha=0.5)
plt.show()
Kent Tower Hamlets
But what about Devon, looking at the list below, didn't we get it all day on the second day?
for visit in r['visits']:
if visit['day'] == 2:
print visit['placeName']
Brixton Yealmpton Modbury Kingsbridge West Charleton Chillington Torcross Stoke Fleming Dartmouth Totnes Paignton Torquay Teignmouth Exeter
We can iterate over the point values and find out which county they are in:
for visit in r['visits']:
if visit['day'] == 2:
for county in uk_counties:
point = visit['Longitude'], visit['Latitude']
if county.geometry.contains(sgeom.Point(point)):
print '{!s:15}\t{}\t{}'.format(visit['placeName'], point, county.attributes['name'])
Brixton (-4.03353, 50.35264) Devon Yealmpton (-4.00218, 50.34773) Devon Modbury (-3.89024, 50.34958) Devon Kingsbridge (-3.77801, 50.28586) Devon West Charleton (-3.74834, 50.26932) Devon Chillington (-3.69358, 50.27124) Devon Torcross (-3.65392, 50.26568) Devon Dartmouth (-3.58328, 50.35206) Devon Totnes (-3.68548, 50.43223) Devon Paignton (-3.56327, 50.4377) Torbay Torquay (-3.52521, 50.472) Torbay Exeter (-3.52429, 50.72216) Devon