I'm attempting to convert latitude/longitude coordinates to X/Y coordinates on a 52,000 sq. ft. grid according to a simple Mercator projection. I used the following Python code to accomplish this, adapted from this post.
from math import pi, log, tan
latitude = 42.374994
longitude = -71.075891
mapWidth = 52000.000000
mapHeight = 52000.000000
x = (longitude+180)*(mapWidth/360)
print "x", x
latRad = latitude*pi/180
print "latRad", latRad
mercN = log(tan((pi/4)+(latRad/2)))
print "mercN", mercN
y = (mapHeight/2)-(mapWidth*mercN/(2*pi))
print "y",y
Using the examples of 30th Street Station and 69th Street Terminal in Philadelphia, I get the following results:
30th Street Station
#39.955641,-75.181178
x 15140.4965111
latRad 0.697357490195
mercN 0.761899319942
y 19694.4777688
69th Street Terminal
#39.962156,-75.258626
x 15129.3095778
latRad 0.697471198396
mercN 0.7620476662
y 19693.2500467
If I calculate the distances between these X/Y coordinates using the Pythagorean theorem, I get a value of
sqrt((15129.30958-15140.49651)^2 + (19693.25005-19694.47777)^2)
= sqrt(125.1474774)
= 11.25410054
I am assuming that in this case the result would take the same units as my mapWidth
and mapHeight
, in feet in this case. But, a value of 11.25 feet is not reasonable. When measured on Google Maps, these two sites are roughly 4 miles apart. What am I missing?