In [38]:
import pandas as pd
import sqlite3
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
from scipy import *
%matplotlib inline
In [2]:
def read_sqllite(sqllite_db, tablename):
""" To read data from sqllite table.
---------------
Parameters
sqllite_db - String. Name of sqllite database
tablename - String. Name of table in database to extract data from
--------------
Returns
df - Dataframe. Data frame of data extracted from database table.
"""
con=sqlite3.connect(sqllite_db)
df=pd.read_sql_query("select * from " + tablename + ";", con)
return df
train_times = read_sqllite(sqllite, 'TIMES')
station_locations = read_sqllite(sqllite, 'LOCATIONS')
In [ ]:
# To extract all stations within an area
d2=d[(d['Easting']>123374) & (d['Easting'] < 340361) & (d['Northing']>6974) & (d['Northing']<152168)]
In [7]:
map2 = Basemap(projection='merc',
resolution = 'l', area_thresh = 0.1,
llcrnrlon=-5.8535212, llcrnrlat=49.898888,
urcrnrlon=-2.856179, urcrnrlat=51.265534)
map2.drawmapboundary(fill_color='aqua')
map2.fillcontinents(color='green', lake_color='aqua')
parallels = np.arange(0.,81,10.)
map2.drawparallels(parallels)
meridians = np.arange(10.,351.,20.)
map2.drawmeridians(meridians
Out[7]:
In [132]:
stations=pd.read_csv('stations.csv')
stations['Easting']=stations['Easting']-80
#stations['Northing']=stations['Northing']-80
In [ ]:
In [133]:
map2 = Basemap(projection='merc',
resolution = 'l', area_thresh = 0.1,
llcrnrlon=-5.8535212, llcrnrlat=49.898888,
urcrnrlon=-2.856179, urcrnrlat=51.265534)
map2.drawmapboundary(fill_color='aqua')
map2.fillcontinents(color='green', lake_color='aqua')
parallels = np.arange(0.,81,10.)
map2.drawparallels(parallels)
meridians = np.arange(10.,351.,20.)
map2.drawmeridians(meridians)
#stations=stations[time==time]
Easting = stations['Easting']
Northing = stations['Northing']
#x,y = map(list(Easting), list(Northing)) # To convert lat/long to projection
#map2.plot(x,y, 'rs', markersize=18)
map2.plot(Easting, Northing, 'rs', markersize=5)
#labels = ['The Needles','Brighton']
#label_style=['rs','bo']
#for label, xpt, ypt, style in zip(labels, x, y, label_style):
# print label, xpt, ypt, style
# map2.plot(xpt,ypt,style, markersize=12)
# plt.text(xpt,ypt, label)
#plt.show
Out[133]:
In [100]:
def OSGB36toWGS84(df):
lats=[]
longs=[]
for station in range(1,len(df), 1):
E = df['Easting'][station]
N = df['Northing'][station]
#E, N are the British national grid coordinates - eastings and northings
a, b = 6377563.396, 6356256.909 #The Airy 180 semi-major and semi-minor axes used for OSGB36 (m)
F0 = 0.9996012717 #scale factor on the central meridian
lat0 = 49*pi/180#Latitude of true origin (radians)
lon0 = -2*pi/180#Longtitude of true origin and central meridian (radians)
N0, E0 = -100000, 400000#Northing & easting of true origin (m)
e2 = 1 - (b*b)/(a*a)#eccentricity squared
n = (a-b)/(a+b)
#Initialise the iterative variables
lat,M = lat0, 0
while N-N0-M >= 0.00001: #Accurate to 0.01mm
lat = (N-N0-M)/(a*F0) + lat;
M1 = (1 + n + (5/4)*n**2 + (5/4)*n**3) * (lat-lat0)
M2 = (3*n + 3*n**2 + (21/8)*n**3) * sin(lat-lat0) * cos(lat+lat0)
M3 = ((15/8)*n**2 + (15/8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
M4 = (35/24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))
#meridional arc
M = b * F0 * (M1 - M2 + M3 - M4)
#transverse radius of curvature
nu = a*F0/sqrt(1-e2*sin(lat)**2)
#meridional radius of curvature
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
eta2 = nu/rho-1
secLat = 1./cos(lat)
VII = tan(lat)/(2*rho*nu)
VIII = tan(lat)/(24*rho*nu**3)*(5+3*tan(lat)**2+eta2-9*tan(lat)**2*eta2)
IX = tan(lat)/(720*rho*nu**5)*(61+90*tan(lat)**2+45*tan(lat)**4)
X = secLat/nu
XI = secLat/(6*nu**3)*(nu/rho+2*tan(lat)**2)
XII = secLat/(120*nu**5)*(5+28*tan(lat)**2+24*tan(lat)**4)
XIIA = secLat/(5040*nu**7)*(61+662*tan(lat)**2+1320*tan(lat)**4+720*tan(lat)**6)
dE = E-E0
lat = lat - VII*dE**2 + VIII*dE**4 - IX*dE**6
lon = lon0 + X*dE - XI*dE**3 + XII*dE**5 - XIIA*dE**7
#Convert to degrees
lat = lat*180/pi
lon = lon*180/pi
#print lat
#print long
lats.append(lat)
longs.append(long)
return lat, long
In [72]:
#stations['Latitude'] = [OSGB36toWGS84(e,n) for e, n in stations[['Easting', 'Northing']]]
#stations['Longitude'] = OSGB36toWGS84(stations[['Easting', 'Northing']])[1]
#[OSGB36toWGS84(e,n) for e, n in stations[['Easting', 'Northing']]]
OSGB36toWGS84(stations[['Easting', 'Northing']])
#stations[['Easting', 'Northing']]
Out[72]:
In [33]:
len(stations)
Out[33]:
In [104]:
from scipy import *
import csv
def OSGB36toWGS84(E,N):
#E, N are the British national grid coordinates - eastings and northings
a, b = 6377563.396, 6356256.909 #The Airy 180 semi-major and semi-minor axes used for OSGB36 (m)
F0 = 0.9996012717 #scale factor on the central meridian
lat0 = 49*pi/180#Latitude of true origin (radians)
lon0 = -2*pi/180#Longtitude of true origin and central meridian (radians)
N0, E0 = -100000, 400000#Northing & easting of true origin (m)
e2 = 1 - (b*b)/(a*a)#eccentricity squared
n = (a-b)/(a+b)
#Initialise the iterative variables
lat,M = lat0, 0
while N-N0-M >= 0.00001: #Accurate to 0.01mm
lat = (N-N0-M)/(a*F0) + lat;
M1 = (1 + n + (5/4)*n**2 + (5/4)*n**3) * (lat-lat0)
M2 = (3*n + 3*n**2 + (21/8)*n**3) * sin(lat-lat0) * cos(lat+lat0)
M3 = ((15/8)*n**2 + (15/8)*n**3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
M4 = (35/24)*n**3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))
#meridional arc
M = b * F0 * (M1 - M2 + M3 - M4)
#transverse radius of curvature
nu = a*F0/sqrt(1-e2*sin(lat)**2)
#meridional radius of curvature
rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)
eta2 = nu/rho-1
secLat = 1./cos(lat)
VII = tan(lat)/(2*rho*nu)
VIII = tan(lat)/(24*rho*nu**3)*(5+3*tan(lat)**2+eta2-9*tan(lat)**2*eta2)
IX = tan(lat)/(720*rho*nu**5)*(61+90*tan(lat)**2+45*tan(lat)**4)
X = secLat/nu
XI = secLat/(6*nu**3)*(nu/rho+2*tan(lat)**2)
XII = secLat/(120*nu**5)*(5+28*tan(lat)**2+24*tan(lat)**4)
XIIA = secLat/(5040*nu**7)*(61+662*tan(lat)**2+1320*tan(lat)**4+720*tan(lat)**6)
dE = E-E0
lat = lat - VII*dE**2 + VIII*dE**4 - IX*dE**6
lon = lon0 + X*dE - XI*dE**3 + XII*dE**5 - XIIA*dE**7
#Convert to degrees
lat = (lat*180)/pi
lon = (lon*180)/pi
return lat, lon
In [107]:
OSGB36toWGS84(340361.0,152168.0)
#OSGB36toWGS84(34.0361, 15.2168)
# - E = 340361
# - N = 152168
Out[107]:
In [101]:
#stations[['Latitude', 'Longitude']] = stations[['Easting', 'Northing']].apply(OSGB36toWGS84())
stations['Latitude'] = stations[['Easting', 'Northing']].apply(OSGB36toWGS84)[0]
#stations['Longitude'] = stations[['Easting', 'Northing']].apply(OSGB36toWGS84)[1]
In [95]:
stations['Latitude'] = stations.apply(lambda row: OSGB36toWGS84(row['Easting'], row['Northing'])[0], axis=1)
In [99]:
stations
Out[99]:
In [ ]: