Jtdylktuy

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:

.

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.

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=Rln(tanphi + 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.

I know this is an old question, but in case you're still interested there is a solution that doesn't involve counting up $cos(phi)$: you can use a binary search to find the latitude of a particular box.

Here's some (python-like) pseudo-code for it:

latitude = 0 degrees

delta = 45 degrees

while convert_latitude_to_box_number(latitude) != desired_box_number:

    if desired_box_number > convert_latitude_to_box_number(latitude):

        latitude += delta

    else:

        latitude -= delta

    delta = delta / 2

//end while (also, latitude now corresponds to desired box)

While this solution is O(log(n)) time and may seem better than the O(n) time $cos(phi)$ solution, you should weigh the cost of computing $R ln(cos(phi)+sec(phi))$ for each iteration.

Since n is bound by 90 degrees in the case of a map, you might be better off just sticking with your $cos(phi)$ solution, but I figured this was worth posting to present another approach.