# -*- coding: utf-8 -*-
"""
Created on Mon Aug 12 13:10:17 2019

@author: jmower
"""

import numpy as np
# Project 74 deg W longitude, 42 deg N latitude to x, y using Mercator projection (for sphere)
R = 6378206.4           # Clarke 1866 equatorial radius
lamb0 = 0               # central meridian for projection
lamb = np.radians(-74)  # user's longitude at 74 degrees W 
phi = np.radians(42)    # user's latitude at 42 degrees N
x = R * (lamb - lamb0)  # compute x
y = R * np.log(np.tan(np.pi/4 + phi/2)) # compute y
print('Forward Projecting with Mercator (for sphere)')
print(f'For longitude {np.degrees(lamb)} and latitude {np.degrees(phi)}, x is {x} meters and y is {y} meters')
print(f'For longitude {lamb} and latitude {phi}, x is {x} meters and y is {y} meters')

# Inverse project x, y above to see if we end up with the original longitude and latitude
print('Inverse Projecting from Mercator coordinates to longitude, latitude')
lambFromInverseProj = x / R + lamb0
phiFromInverseProj = np.pi/2 - 2 * np.arctan(np.power(np.e, -y/R))
print (f'for x = {x}, y = {y}, lamb = {np.degrees(lambFromInverseProj)}, phi = {np.degrees(phiFromInverseProj)}')