mercator projections

[jjguy]

If input points are lat/long and cover a reasonably large distance (10s of miles) the resulting points will be slightly off due to the lack of mercator projections.

The translation is simple:

from math import degrees, radians, asinh, tan
oldLatRads = radians(lat)
newLatRads = asinh(tan(oldLatRadians))
newLat = degrees(newLatRads)

This can't go in saveKML(), since that operation is post-image translation. Simplest impl is a new function provided by the lib to translate a set of points before calling heatmap(). This will add one pass in native python, decreasing perf for large sets of lat/long.

[jjguy]

Heh, not that simple. I implemented below, the results are unexpected.

def mercator(self, pts):
    """ 
    Adjust the latitude in a list of coordinate degrees for accurate map layout.

    pts -> an iterable list of lat/long tuples: [(lat,long), (lat,long), ...]

    Returns a new copy of the list.
    """
    for ndx, (lat, lng) in enumerate(pts):
        # mercator projection is asin(tan(lat)) in radians
        newlat = math.degrees(math.asinh(math.tan(math.radians(lat))))
        pts[ndx] = (newlat, lng)

    return pts 

[kwauchope]

You will probably want to use something like https://pypi.python.org/pypi/pyproj/.

I'm also not convinced the translate function is overly accurate over large distances either. I'll have to do some reading.

[kwauchope]

Got a new branch where I'm testing the projections. Requires pyproj. I don't think Heatmap.c will need to be changed as long as everything is done in Cylindrical linear projections or simple x,y so the linear translation still works. Can input whatever projection you like to heatmap.py as long as it's supported by pyproj. Default is no projection at all to be backwards compatible.

[jjguy]

@kwauchope sorry for not giving this more attention recently; balancing a lot in life at the moment.

send a pull request on, I'll merge it in and cut a new point release with support. my only request beyond what you've provided is to make the projection support conditional, based on the pyproj availability. i.e.,

• no projection, no pyproj import
• projection enabled, try import
  -- if success, use pyproj
  -- if not success, print a warning and fall back to direct plotting

thanks!

[kwauchope]

I've added the option for pyproj, wrote the warning as stdout if the input projection is set (if there is no input projection it is ignored anyway). It doesn't work particularly well with the test output though if pyproj is not installed. Not sure of a better way to show the warning though.

Perhaps it can be an outright exception anyway, if a user is expecting the tool to do projections by giving it at least a source projection then shouldn't it just quit and go home?

[kwauchope]

I hate to add more complexity but especially with the allowance of different projections, not all pixels are created as equals so dotsize may not be a reliable way to present the data. Say with cylindrical projections the pixels are worth more metres towards the equator then the poles so it would give higher density to values around the equator compared to the poles. Only a noticeable issue for larger/global datasets but it is still not ideal.

A more reliable measure of size for geospatial projections would be degrees or metres. This means that each calculation in heatmap.c would have different step sizes. In addition to the projected x,y and optional weight there would be a separate pixels measurement for each point. Not hard to fix though adds extra overhead. Other things to work on at the moment but it is something I will look at fixing in the future for a really robust solution,

[kwauchope]

@jjguy I was wondering your thoughts on disabling KML output if pyproj isn't installed. Say the input is WGS84 (the format KML needs for the overlay bounding box and the most probable input format) the linear interpolation in the image will be wrong. If the input is World Equidistant Cylindrical the linear interpolation will be correct but bounding box will be out.