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
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?