[FEATURE] add built in coordinate conversions

[comperem]

Hello,

The GIS and mapping capability of Bokeh is good and growing but lacks ability to convert units between (lat,lon) pairs and Web Mercator Easting and Northing coordinate pairs.

pyproj is a simple, focused and mature coordinate conversion utility for converting between recognized standard coordinate frames.

Adding pyproj or a similar capability would allow native Bokeh code and examples to perform coordinate conversions as a normal part of GIS mapping and plotting. For example, this code and associated documentation example could be included:

I propose adding pyproj. How can this get considered for the next version update?

Alternative codes and methods:

1. OpenStreetMap Mercator wiki with conversions written in multiple programming languages
2. Web Mercator projection wiki has formulas and good explanations

References:

1. discourse link located here.
2. pull request with working example except for pyproj dependency.

[bryevdv]

cc @bokeh/dev for for input

I'll note that historically we have tried to keep Bokeh's dependency list very slim. That said, the current situation where folks have to convert their own web mercator is indeed very bad UX for users. @comperem can you confirm that pyproj is available on conda and PyPI?

[jbednar]

For what it's worth, Datashader also needed the specific projection from lon,lat to web Mercator for our examples, and to avoid having a dependency on pyproj we just added that specific function in a numpy/pandas compatible form:

import numpy as np 
def lnglat_to_meters(longitude, latitude):
    """
    Projects the given (longitude, latitude) values into Web Mercator
    coordinates (meters East of Greenwich and meters North of the Equator).

    Longitude and latitude can be provided as scalars, Pandas columns,
    or Numpy arrays, and will be returned in the same form.  Lists
    or tuples will be converted to Numpy arrays.

    Examples:
       easting, northing = lnglat_to_meters(-40.71,74)

       easting, northing = lnglat_to_meters(np.array([-74]),np.array([40.71]))

       df=pandas.DataFrame(dict(longitude=np.array([-74]),latitude=np.array([40.71])))
       df.loc[:, 'longitude'], df.loc[:, 'latitude'] = lnglat_to_meters(df.longitude,df.latitude)
    """
    if isinstance(longitude, (list, tuple)):
        longitude = np.array(longitude)
    if isinstance(latitude, (list, tuple)):
        latitude = np.array(latitude)

    origin_shift = np.pi * 6378137
    easting = longitude * origin_shift / 180.0
    northing = np.log(np.tan((90 + latitude) * np.pi / 360.0)) * origin_shift / np.pi
    return (easting, northing)

[comperem]

@bryevdv , yes it appears Pyproj is available via both PyPi and Conda.

The Pyproj installation instructions describe both here: http://pyproj4.github.io/pyproj/stable/installation.html

The PyPi link is here: https://pypi.org/project/pyproj/

The Conda source is a little less direct because it's on the Conda-forge channel, but their install line specifies this so it's still a simple install: conda install -c conda-forge pyproj

It appears pretty active and well maintained.

[philippjfr]

I see you've closed this, but I'll just chime in to say that pyproj is not a trivial dependency and I would certainly recommend against requiring it by default particularly if the aim is merely to provide lat/lon -> Mercator projections. As @jbednar points out the same can be achieved with a fairly trivial function.

[bryevdv]

Right conda-forge only is a big negative. I'll re-open this for adding the function above. Where is the best place to put it?

[comperem]

@jbednar , this sounds straightforward, and it seems you've gotten a good simple solution, which is good.

But, I rarely do a Group oriented project without needing UTM as well.

[comperem]

this seems like a short-sighted evaluation. Web Mercator is not the only need because:

(a) there are more unit conversions commonly used than (lat,lon) to Mercator, and

(b) there is substantial detail and technique that ensure accuracy and standards compliance for the different world sphere or ellipsoid assumptions and datums associated with WGS-84 and others.

World geodetic systems are a complex mix of international standards and industry practice. The proj library, presumably, stays aware and current on these issues. pyproj will provide that expertise over time with clear documentation.

Web Mercator is prevalent in all the main map sites like Google and Bing and others. So it's critical for display but not convenient for distance or area calculations. It's a global projection and not an orthogonal Cartesian coordinate frame so distances are not straightforward.

UTM is quite useful for achieving units of meters from a (lat,lon) pair. The globe is divided into 60 UTM zones with each assumed to be an orthogonal XYZ coordinate frame for that whole slice. You can calculate distance and areas within each zone just like you'd expect: the square root of the sum of the squares of the deltas in X, Y, and Z.

So we see, (lat,lon), Web Mercator and UTM are all necessary for mapping and associated maths.

Having ability to move easily between these is kind of mandatory.

Inclusion in Bokeh is a different story. I use pyproj and utm in most mapping projects.

If there's a way to use these in Bokeh examples without tying them as Bokeh native dependencies that might be easier. I just don't know how to accomplish that.

Somehow making these coordinate conversion tools available with Bokeh seems to make a lot of sense and would facilitate more GIS use with Bokeh.

[comperem]

that one function may be good enough to make good examples to illustrate what's possible with Bokeh maps.

Anyone who wants more can include more in their own projects.

That one function may be the simplest approach from a software maintenance standpoint.

I haven't studied that particular function but @jbednar says it's provided good results.

[bryevdv]

Well, it's late, I was a bit circumspect. I guess a more thorough exposition is: we should add our own functions, whatever those might be. Bokeh really always needs to be able to be installed from defaults, i.e. with just conda install bokeh So a hard dependency on something that is only on conda-forge is a non-starter, unfortunately.

[comperem]

I get it. Any adjustment cannot cause install problems.

[comperem]

UTM is something needed for additional work beyond the basic demo in this initial post.

Lat, lon to Web Mercator is the only conversion needed to get that demo working.

[bryevdv]

@comperem Do you think you would you be able to add the minimal conversion functions to make the demo better (and added to the docs) by the end of next week? Would love to see this in 2.1 and that is our code freeze cutoff. And then we can continue to discuss additional transformations (either through vendoring, or less onerous dependencies) here.

Offhand I suppose a new bokeh.geo module is probably the best place for these to go.

[comperem]

Bryan,

Yes, I'll try to put something together. The function posted earlier may very well be sufficient. I found another function that is not dependent upon Numpy at OpenStreetMaps. But a Numpy option may be helpful for some circumstances. I'll try both to compare both and provide simple, clear, minimal documentation.

By the end of this upcoming week.

Marc

On Thu, May 21, 2020, 11:45 AM Bryan Van de Ven notifications@github.com wrote:

@comperem https://github.com/comperem Do you think you would you be able to add the minimal conversion functions to make the demo better (and added to the docs) by the end of next week? Would love to see this in 2.1 and that is our code freeze cutoff. And then we can continue to discuss additional transformations (either through vendoring, or less onerous dependencies) here.

Offhand I suppose a new bokeh.geo module is probably the best place for these to go.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/bokeh/bokeh/issues/10009#issuecomment-632167476, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG7MPYUX2WE2NL5UF27XI4TRSVECPANCNFSM4M6IT2WQ .

[bryevdv]

@comperem FYI Numpy is a hard requirement of Bokeh these days so it is fine to use

[comperem]

@bryevdv, The code listed above from @jbednar looks pretty good and compares identically to the comparable coordinate conversion using pyproj(). The coordinate conversion changes (lat/lon) in decimal degrees to Web Mercator Easting and Northing values.

The only correction I see necessary is a typo in the example conversion documented inside his lnglat_to_meters() function. They seem to be listed backwards and negative on both. So switching the first Examples: to this should make the documentation consistent across all three examples provided:

easting, northing = lnglat_to_meters(-74,40.71)

Otherwise, his function is helpful. It would make mapping functions work easily in native Bokeh.

[jbednar]

You're right; the examples would be clearer if they all used the same point! The updated docstring is below, and I've fixed it on Datashader as well.

import numpy as np 
def lnglat_to_meters(longitude, latitude):
    """
    Projects the given (longitude, latitude) values into Web Mercator
    coordinates (meters East of Greenwich and meters North of the Equator).

    Longitude and latitude can be provided as scalars, Pandas columns,
    or Numpy arrays, and will be returned in the same form.  Lists
    or tuples will be converted to Numpy arrays.

    Examples:
       easting, northing = lnglat_to_meters(-74,40.71)

       easting, northing = lnglat_to_meters(np.array([-74]),np.array([40.71]))

       df=pandas.DataFrame(dict(longitude=np.array([-74]),latitude=np.array([40.71])))
       df.loc[:, 'longitude'], df.loc[:, 'latitude'] = lnglat_to_meters(df.longitude,df.latitude)
    """
    if isinstance(longitude, (list, tuple)):
        longitude = np.array(longitude)
    if isinstance(latitude, (list, tuple)):
        latitude = np.array(latitude)

    origin_shift = np.pi * 6378137
    easting = longitude * origin_shift / 180.0
    northing = np.log(np.tan((90 + latitude) * np.pi / 360.0)) * origin_shift / np.pi
    return (easting, northing)

[comperem]

@jbednar, it seems the native python function you've shown here is about 4 orders of magnitude faster than the Web Mercator coordinate conversion with pyproj. My code had a bottleneck and was limited by the pyproj conversion. Timing tests between the two showed approximately 14,000 times speed improvement so I'll be using your function from now on. Thanks

[bryevdv]

FWIW I am definitely +1 on adding simple functions to convert from lat/lng to meters and vice versa. Someone just needs to make some specific proposals about where, etc. This would be a great self-contained task for a new contributor.

[jbednar]

14,000 times speed improvement

Whoa! I'm glad it was useful. The code is Pythonic but the heavy lifting is all done in C by NumPy, so Python doesn't deserve all the credit. :-)

[comperem]

It was surprising. I had a udp buffer overflow and tracked down the delay to pyproj. I have no idea what they're doing to take so long inside that library. Python, even if it's an interface, and C are such a powerful combination.

On Fri, Aug 7, 2020 at 10:09 PM James A. Bednar notifications@github.com wrote:

14,000 times speed improvement

Whoa! I'm glad it was useful. The code is Pythonic but the heavy lifting is all done in C by NumPy, so Python doesn't deserve all the credit. :-)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/bokeh/bokeh/issues/10009#issuecomment-670808472, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG7MPYS7NSRFJVW7DTUNKFLR7SXWLANCNFSM4M6IT2WQ .

[comperem]

@bryevdv ,

Did the simple GIS coordinate conversion submitted by @jbednar [immediately above] (https://github.com/bokeh/bokeh/issues/10009#issuecomment-636555037) ever get added?

It looks like this task was dropped, or never initiated. If I recall correctly, this unit conversion code was the only remaining detail needed to get the mapping tile demo integrated into the Bokeh code base.

What's the next step to completion?

This thread is already a feature request. How about putting that function in a new file in a new folder called GIS_utils or something similar?

Nobody knows how to proceed except seasoned maintainers.

Marc