[Math] Simpler mathematic formula to find latitude coordinate mapping to lines “equally sized” on mercator projection

geodesytrigonometry

I'm implementing a map visualization atop a mercator projected map (e.g google maps) where each circle appears to be the same size on the map:

enter image description here.

At the equator, a circle plotted with a one degree latitude x one degree longitude looks like a circle, but as you go further away from the equator, the height of a one degree appears taller on the map, so we scale it back. Thus the "equally sized" circles are in fact squatter and squatter dimensions.

We're taking a number of lat/lon points and mapping them to the appropriate circle, and then plotting the circles with darker colors indicating a higher density. So this finally leads to may question: is there a simple formula to calculate the correct latitude offset for a given latitude point? For longitude, it's simply the floor, but for latitude, this won't work.

I'm currently finding the appropriate 'bin' or floor for the latitude by counting up cos(angle) size chunks from the equator; at zero, one degree works and as you go north, the delta decreases. Here's sample python code:

def makelatbins():
    """
    make steps of latitude that will appear equal hieght on the
    mercator projecion based google map; if we did equal size degrees
    the height of each circle gets taller as you go further north.
    """
    import math
    latbins = []
    lat = 0
    while lat < 84:
        latbins.append(lat)
        lat += math.cos(math.radians(lat))
    return latbins

latbins = makelatbins()

def latbin(lat):
    import bisect
    index = bisect.bisect_left(latbins, lat)
    return latbins[index]

I have this nagging sense that there should be a simple formula to compute the appropriate floor function for a latitude point, but can't figure it out.

Best Answer

On a Mercator projection the verticals are lines of constant longitude with equal spacing. The horizontals are lines of constant latitude, with $x=R(\lambda-\lambda_0)$, where $R$ is the radius of the earth and $\lambda_0$ is the central longitude of the map. $y=R\ln(\tan\phi + \sec \phi)$ so a small step in $y$ corresponds to $R \sec \phi$ in distance. So you should scale your verticals by $1/\sec \phi$ to keep the circles round on the map.

Related Question