Convert lat/long coordinates into an ascii string. Useful for database storage outside of geodata formats.

encodeLongLatCoords.py

# Objective: A system capable of encoding geographic points, and polylines in ASCII characters
# 			 Obviously Google and geohash have similar algorithms, but with limited precision
#
# A 24-bit integer can store binary integer values up to 16,777,215
# Using PEM format, 24 bits can be represented as 4 base-64 characters with the alphabet:
# ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/
# A polyline can be represented by a string of characters,
# where each 8 characters is a longitude/latitude location
# For maximum flexibility, a polyline or point string is prefixed with three characters
# The first two characters specifies the origin - the lower left corner of a MGRS GZD
# The MGRS longitude value (zero to sixty) will be replaced with its base-64 equivalent
# The third character is a code of precision:
#		d (22): double precision to about 0.0000002 arc degrees (about 1 inch)
#		f (21): full precision (default) to about 0.0000005 arc degrees (about 5 cm)
#		h (20): half precision to 0.000001 arc degree (about 11 cm)
#		q (19): quarter precision to 0.000002 arc degrees (about 22 cm)
#		e (18): eighth precision to 0.000004 arc degrees (about half meter)
#		x (17): sixteenth precision to 0.000008 arc degrees (about 1 meter)
#		t (16): 1/32nd precision to 0.000016 arc degrees (about 2 meters)
#		s (15): 1/64th precision to 0.000032 arc degrees (about 4 meters)
# At resolution "s", a polyline could be encoded that encircles the globe
# Google polyline encoding is precise to 0.00001, which is roughtly equivalent to 'x' (1 m).
# Precision "d" is not very useful because it covers only the lower left quarter of an MGRS GZD
# Future update could refine the MGRS system to allow for double precision
#
# Example: Encoded point = 'JShRy45aPAZ'
# First three characters 'JSh' indicates MGRS GZD of 10S, with half precision (J=10, S=S, h=h)
# The first coordinate is longitude 'Ry45' by latitude 'aPAZ', which are 24-bit integer values
# encoded in base-64. 'Ry45' is three bytes: 71, 46, 57 or interpreted as one 24-bit integer
# 'Ry45' equals 4,664,889, which represents 4664889 / 2**20 = 4.448785 and 'aPAZ' equals 6.558618.
# The power of two ('h'=20) is used from the precision to convert the integer to a decimal number
# These coordinates when added to the origin of MGRS GZD 10S (-126,32) gives -121.551215, 38.558618
# The exact same point could also be coded JSfjlxy0eAz, with more precision.

import base64
from bisect import bisect

def decodeMGRS(code):
	b64st = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz01234567'
	mgrsLatCodes = 'CDEFGHJKLMNPQRSTUVWX'
	mgrsLat = {c:l for c,l in zip(mgrsLatCodes,range(-80,80,8))}
	mgrsLong = {c:l for c,l in zip(b64st,range(-180,180,6))}
	return (mgrsLong[code[0]], mgrsLat[code[1]])

def encodeMGRS(longLat):
	b64st = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz01234567'
	mgrsLatCodes = 'CDEFGHJKLMNPQRSTUVWX'
	mgrsLat = {l:c for c,l in zip(mgrsLatCodes,range(-80,80,8))}
	mgrsLong = {l:c for c,l in zip(b64st,range(-180,180,6))}
	return mgrsLong[longLat[0]] + mgrsLat[longLat[1]]

def findMGRSorigin(longLat):
	LONG, LAT = longLat
	MGRSlongs = [x for x in range(-180,180,6)]
	MGRSlats = [y for y in range(-80,80,8)]
	x = MGRSlongs[bisect(MGRSlongs,LONG)-1]
	y = MGRSlats[bisect(MGRSlats,LAT)-1]
	return encodeMGRS((x, y))
	
def encode24bitInt(val):
	b3 = int(val % 256)
	b2 = int(((val - b3)/256) % 256)
	b1 = int((val - b2*256 - b3)/256/256)
	return base64.b64encode(bytearray([b1,b2,b3])).decode()

def decode24bitInt(str):
	b1, b2, b3 = base64.b64decode(str)
	return b3 + b2*256 + b1*256*256

def roundInt(decimal):
	return int(round(decimal,0))

enc = encode24bitInt(15589981)
print (enc + ' is equal to %i' % decode24bitInt(enc))

def encodeLongLat(coordinates, origin='auto', precision=21):
	preCode = '012345678901234stxeqhfd'
	long, lat = coordinates
	if origin == 'auto':
		origin = findMGRSorigin(coordinates)
	elif type(origin) == tuple or type(origin) == list:
		oLong, oLat = origin
		origin = encodeMGRS(origin)
	if type(origin) == str:
		oLong, oLat = decodeMGRS(origin)
	long, lat = (long-oLong)*2**precision, (lat-oLat)*2**precision
	return (origin +preCode[precision] + 
			encode24bitInt(roundInt(long)) + 
			encode24bitInt(roundInt(lat)))

def decodeLongLat(coordString):
	preCode = '012345678901234stxeqhfd'
	origin = decodeMGRS(coordString[:2])
	precision = preCode.find(coordString[2])
	long = coordString[3:7]
	lat = coordString[7:]
	return (float(decode24bitInt(long))/2**precision+origin[0], 
		float(decode24bitInt(lat))/2**precision+origin[1])

enc2 = encodeLongLat((-122.5512156, 36.5586184))
print (enc2 + ' is equal to %f, %f' % decodeLongLat(enc2))

def extent(coordList):
	xMin, yMin = 180, 90
	xMax, yMax = -180, -90
	for coord in coordList:
		xMin = min(xMin,coord[0])
		xMax = max(xMax,coord[0])
		yMin = min(yMin,coord[1])
		yMax = max(yMax,coord[1])
	return ((xMin,yMin),(xMax,yMax))

def maxPrecision(origin, maxXY):
	xLength = maxXY[0] - origin[0]
	yLength = maxXY[1] - origin[1]
	maxLength = max(xLength, yLength)
	cutoffs = [8,16,32,64,128,256,512]
	precisions = [21,20,19,18,17,16,15]
	return precisions[bisect(cutoffs,maxLength)]

def encodePolyline(coordList, origin='auto', precision='auto'):
	preCode = '012345678901234stxeqhfd'
	if origin == 'auto' or precision == 'auto':
		minXY, maxXY = extent(coordList)
	if origin == 'auto':
		origin = findMGRSorigin(minXY)
	if precision == 'auto':
		precision = maxPrecision(decodeMGRS(origin),maxXY)
		print (precision)
	coordString = "".join([s[3:] for s in 
						[encodeLongLat(coords,origin,precision) for 
						coords in coordList]])
	return origin + preCode[precision] + coordString

def decodePolyline(lineString):
	coordList = []
	numCoords = int((len(lineString)-3)/8)
	origPrecis = lineString[:3]
	for c in range(numCoords):
		coordList.append(decodeLongLat(origPrecis+lineString[c*8+3:c*8+11]))
	return coordList

pl = [(-121.551215,38.558618),(-121.569987,38.670099),(-121.579908,41.789957)]
enc3 = encodePolyline(pl)
print (enc3 + ' translates to the line: ')
print (decodePolyline(enc3))