Transforming geodesic latitude, longitude coordinates to geocentric X, Y, Z coordinates and back again.

ellipsoid.py

'''
Define an ellipsoid by its semi major axis (a) and flattening (f).
The definition of the ellipsoid is then saved in the instance object and the functions can be called with just the coordinates to be transformed.
'''

from math import *

class Ellipsoid:
	def __init__(self, a, f):
		self.a = a
		self.f = 1/f
		self.eSq = self.geteSquared(self.f)

	def getX(self, lat, lon, h):
		return (self.getNprime(lat) + h) * cos(radians(lat)) * cos(radians(lon))

	def getY(self, lat, lon, h):
		return (self.getNprime(lat) + h) * cos(radians(lat)) * sin(radians(lon))

	def getZ(self, lat, h):
		return (self.getNprime(lat) * (1 - self.eSq) + h) * sin(radians(lat))

	def getGeocentric(self, lat, lon, h):
		return [self.getX(lat, lon, h), self.getY(lat, lon, h), self.getZ(lat, h)]

	def getNprime(self, lat):
		return self.a / sqrt(1 - self.eSq * pow(sin(radians(lat)), 2))

	def geteSquared(self, f):
		return f * (2 - f)

	def getLatLon(self, geocentric):
		lon = degrees(atan(geocentric[1] / geocentric[0]))
		p = sqrt((geocentric[0]**2) + (geocentric[1]**2))
		theta = degrees(atan((geocentric[2]) / (p * sqrt(1 - self.eSq))))
		# prep for lat
		inner_division = (self.a * self.eSq) / (sqrt(1 - self.eSq))
		above_division = geocentric[2] + inner_division * (sin(radians(theta)) ** 3)
		#below division
		below_division = p - self.a * self.eSq * (cos(radians(theta)) ** 3)
		lat = degrees(atan(above_division / below_division))
		h = ((p) / (cos(radians(lat)))) - self.getNprime(lat)
		return [lat, lon, h]

geocentric.py

'''
Johnnie Hard, 2017
johnnie.hard@gmail.com
Coursework for GE7030, Positioning, Map Projections and Digital Photogrammetry, 7.5 hp
Master's Programme in Geomatics with Remote Sensing and GIS, Stockholm University

Exercise:
From a given lat, lon pair (and height above the ellipsoid), calculate the geocentric XYZ coordinates and then back to lat, lon.

Formulas: http://www.lantmateriet.se/globalassets/kartor-och-geografisk-information/gps-och-matning/geodesi/transformationer/xyz_geodetiska_koord_och_exempel.pdf
'''
from ellipsoid import *

print "\nInitialize grs80 and bessel ellipsoids...\n"
grs80 = Ellipsoid(6378137, 298.257222101)
bessel = Ellipsoid(6377397.155, 299.1528128)

print "Get geocentric X, Y, Z for lat, lon, h coordinates 58, 17, 30...\n"
grs80_geocentric = grs80.getGeocentric(58, 17, 30)
bessel_geocentric = bessel.getGeocentric(58, 17, 30)

print 'grs80: ', 'X', grs80_geocentric[0], 'Y', grs80_geocentric[1], 'Z', grs80_geocentric[2]
print 'bessel: ', 'X', bessel_geocentric[0], 'Y', bessel_geocentric[1], 'Z', bessel_geocentric[2]

print '\nAnd back to lat, lon:\n'

grs80_latlon = grs80.getLatLon(grs80_geocentric)
bessel_latlon = bessel.getLatLon(bessel_geocentric)

print 'grs80: ', grs80_latlon[0], grs80_latlon[1], grs80_latlon[2]
print 'bessel: ', bessel_latlon[0], bessel_latlon[1], bessel_latlon[2]