import csv

from numpy import pi,sin,cos,sqrt,tan,arctan2

def read_csv(fname='data/barnet_data2.csv',header=True):

with open(fname,"r") as f:

reader = csv.reader(f)

if header:

next(reader,None)

data = []

for row in reader:

data.append(row)

return data

def writeListToFile(list,file='data/output.txt'):

with open(file, 'w') as f:

for item in list:

f.write("%f\n" % item)

def write2csv(data,file):

# strdata = str(data)

strdata = ','.join('%s' % str(x) for x in data)

with open(file, 'a') as f:

f.write(strdata+'\n')

# http://hannahfry.co.uk/2012/02/01/converting-latitude-and-longitude-to-british-national-grid/

def WGS84toOSGB36(lat, lon):

#First convert to radians

#These are on the wrong ellipsoid currently: GRS80. (Denoted by _1)

lat_1 = lat*pi/180

lon_1 = lon*pi/180

#Want to convert to the Airy 1830 ellipsoid, which has the following:

a_1, b_1 =6378137.000, 6356752.3141 #The GSR80 semi-major and semi-minor axes used for WGS84(m)

e2_1 = 1- (b_1*b_1)/(a_1*a_1) #The eccentricity of the GRS80 ellipsoid

nu_1 = a_1/sqrt(1-e2_1*sin(lat_1)**2)

#First convert to cartesian from spherical polar coordinates

H = 0 #Third spherical coord.

x_1 = (nu_1 + H)*cos(lat_1)*cos(lon_1)

y_1 = (nu_1+ H)*cos(lat_1)*sin(lon_1)

z_1 = ((1-e2_1)*nu_1 +H)*sin(lat_1)

#Perform Helmut transform (to go between GRS80 (_1) and Airy 1830 (_2))

s = 20.4894*10**-6 #The scale factor -1

tx, ty, tz = -446.448, 125.157, -542.060 #The translations along x,y,z axes respectively

rxs,rys,rzs = -0.1502, -0.2470, -0.8421 #The rotations along x,y,z respectively, in seconds

rx, ry, rz = rxs*pi/(180*3600.), rys*pi/(180*3600.), rzs*pi/(180*3600.) #In radians

x_2 = tx + (1+s)*x_1 + (-rz)*y_1 + (ry)*z_1

y_2 = ty + (rz)*x_1 + (1+s)*y_1 + (-rx)*z_1

z_2 = tz + (-ry)*x_1 + (rx)*y_1 + (1+s)*z_1

#Back to spherical polar coordinates from cartesian

#Need some of the characteristics of the new ellipsoid

a, b = 6377563.396, 6356256.909 #The GSR80 semi-major and semi-minor axes used for WGS84(m)

e2 = 1- (b*b)/(a*a) #The eccentricity of the Airy 1830 ellipsoid

p = sqrt(x_2**2 + y_2**2)

#Lat is obtained by an iterative proceedure:

lat = arctan2(z_2,(p*(1-e2))) #Initial value

latold = 2*pi

while abs(lat - latold)>10**-16:

lat, latold = latold, lat

nu = a/sqrt(1-e2*sin(latold)**2)

lat = arctan2(z_2+e2*nu*sin(latold), p)

#Lon and height are then pretty easy

lon = arctan2(y_2,x_2)

H = p/cos(lat) - nu

#E, N are the British national grid coordinates - eastings and northings

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)

n = (a-b)/(a+b)

#meridional radius of curvature

rho = a*F0*(1-e2)*(1-e2*sin(lat)**2)**(-1.5)

eta2 = nu*F0/rho-1

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)

I = M + N0

II = nu*F0*sin(lat)*cos(lat)/2

III = nu*F0*sin(lat)*cos(lat)**3*(5- tan(lat)**2 + 9*eta2)/24

IIIA = nu*F0*sin(lat)*cos(lat)**5*(61- 58*tan(lat)**2 + tan(lat)**4)/720

IV = nu*F0*cos(lat)

V = nu*F0*cos(lat)**3*(nu/rho - tan(lat)**2)/6

VI = nu*F0*cos(lat)**5*(5 - 18* tan(lat)**2 + tan(lat)**4 + 14*eta2 - 58*eta2*tan(lat)**2)/120

N = I + II*(lon-lon0)**2 + III*(lon- lon0)**4 + IIIA*(lon-lon0)**6

E = E0 + IV*(lon-lon0) + V*(lon- lon0)**3 + VI*(lon- lon0)**5

#Job's a good'n.

return E,N