Convert lat/lon to Mercator x/y, then measure distance

Question

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?

Answer 1

You have projected the whole globe to a 52000 x 52000 square. Obviously, the unit here cannot be feet as, 52000 feet would be equal to ca. 9.85 miles. But speaking of unit, it cannot really be defined in such length as on different parts of the maps it would mean different "real world" length. As per the grid:

"... on a 52,000 sq. ft. grid ..."

This is false, your grid is 2 704 000 000 sq. ft., as 52000 ft is one side of it.