1
$\begingroup$

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.

2 Answers 2

2

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.

  • 0
    Yes, scaling by 1/secϕ, or cos(ϕ) gives the correct ratio, but in addition to finding the correct height of the circle, I need to know the latitude coordinate that maps to an equal height circles all the way from 0 latitude. That's why right now I'm counting in discreet steps of cos(ϕ) from ϕ = 0. But I'm wondering if there's a formulaic way to find the nearest boundary without counting up from 0.2011-10-27
  • 0
    The equation $y=R\ln(\tan\phi + \sec \phi)$ shows how to do that, as it is the integral of $\sec \phi$. $R$ is not the earth radius, it is just your scale factor. To find what $\phi$ corresponds to a particular $y$, you need to invert this, which you can do numerically.2011-10-27
  • 0
    I think I'm getting closer to understanding this, perhaps a concrete example will help. I'm mapping all of a collection of lat,lons to bins of equal square sizes on the map, and then I will plot circles for each square. For a point at position 78.3, 68.3, I need to find which counter this maps to. I can map it to longitude position 68.0, since each longitude box is equal width. I know the height of the box will be scaled by cos(78.3) = 0.203. But where is the starting point of the box in latitude?2011-10-27
  • 0
    (continued) You are saying that I can use y=R*ln(tan(lat) + sec(lat)) to find the starting point of the box. What is R here? The scaling factor, cos(78.3)? (upvoting for the time you have taken already).2011-10-27
  • 0
    You have defined your horizontal scale to be one degree per box, so if you plot the whole earth you need to have a width of 360 boxes. In this case $R=180/\pi$ to make the $x$ and $y$ scales match as well as possible. In the other axis you need to decide how far north you are going to go. If your map is square, it will go up to 85.05113 deg latitude, corresponding to $y=\pi$ for 180 boxes. Latitude 78.3 plots to $y=2.2783$, or box number 130.5372011-10-27
  • 0
    Ok Great! I can indeed confirm that 78.3 maps into the 130th box from the equator, and it can be calculated via floor((180 / pi) * (ln(tan(78.3) + sec(78.3)))). So the final step is, how do I figure out the bottom coordinate of box number 130 without counting up cos(lat) at a time?2011-10-27
  • 0
    I don't think you can get the exact value without a numeric approach. You would wind up solving $\tan \phi+\sec \phi=\frac{1+\sin \phi}{\cos \phi}=\exp(y\pi/180)$2011-10-27
1

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.