Fix axis tick labels for non-Mercator projections

[suvarchal]

I use the functionality of bokeh tools in my plots comprising coastlines, the default Web-Mercator projection makes my plots not very appealing. So in almost all of my applications, I use the following function to get a rectangular projection/native projection of the data. I wonder if something like this should go inside geoviews.feature when bokeh backend is used or at least part of bokeh examples. This was limiting for me for a while so I guess should be useful for others.

def coastlines(resolution='110m',lon_360=False):
    """ A custom method to plot in cylyndrical equi projection, most useful for
        native projections, geoviews currently supports only Web Mercator in
        bokeh mode.
        
        Other resolutions can be 50m
        
        lon_360 flag specifies if longitudes are from -180 to 180 (default) or 0 to 360
        TODO: if hv.Polygons is used instead of overlay it is way faster but 
              something is wrong there.
    """
    try:
        import cartopy.io.shapereader as shapereader
        from cartopy.io.shapereader import natural_earth
        import shapefile
        filename = natural_earth(resolution=resolution,category='physical',name='coastline')
    
        sf = shapefile.Reader(filename)
        fields = [x[0] for x in sf.fields][1:]
        records = sf.records()
        shps = [s.points for s in sf.shapes()]
        pls=[]
        for shp in shps:
            lons=[lo for lo,_ in shp]
            lats=[la for _,la in shp]
            if lon_360:
                lats_patch1=[lat for lon,lat in zip(lons,lats) if lon<0]
                lons_patch1=[lon+360.0 for lon in lons if lon<0]
                if any(lons_patch1):
                    pls.append(hv.Path((lons_patch1,lats_patch1))(style={'color':'Black'}))    
                lats_patch2=[lat if lon>=0 else None for lon,lat in zip(lons,lats)]
                lons_patch2=[lon if lon>=0 else None for lon in lons]
                if any(lons_patch2):
                    pls.append(hv.Path((lons_patch2,lats_patch2))(style={'color':'Black'}))
            else:
                pls.append(hv.Path((lons,lats))(style={'color':'Black'}))
        return hv.Overlay(pls)
    except Exception as err:
        print('Overlaying Coastlines not available from holoviews because: {0}'.format(err))

#example
coastline=coastlines(lon_360=True) #this may download shape files on first invocation
                                   #if data longitudes are -180 to 180 then lon_360=False

#use it on holoviews elements
img*coastline

not sure if it is more holoviews or geoviews related, but you are a better judge if you think it is useful.

[philippjfr]

@suvarchal GeoViews has supported plotting arbitrary projections with bokeh since the 1.5 release, see http://blog.holoviews.org/release_1.5.html and http://geoviews.org/user_guide/Projections.html. So you can now do something like:

coastline = gv.feature.coastline.options(projection=ccrs.PlateCarree(central_longitude=180), global_extent=True, width=600)

coastline * gv.Points(np.random.randn(100, 2)*50)

So that works now. That said you'll also notice that while the points are in the right space the axes aren't labelled quite right, so if you don't mind I'll repurpose this issue to address that.

[suvarchal]

Cool, I have been using this for months now and meanwhile haven't caught up with recent changes in geoviews (early adapter syndrome!).

[suvarchal]

Sure

So that works now. That said you'll also notice that while the points are in the right space the axes aren't labelled quite right, so if you don't mind I'll repurpose this issue to address that.

[philippjfr]

Not sure we can handle this in a totally general way for all projections but at minimum we should at least try to fix it for PlateCarree and RotatedPole cases.

[suvarchal]

Exactly, the most used cases would be the most useful.

[jbednar]

If you're referring to longitude being offset by 180 degrees, what's the underlying cause of that? It seems strange that it wouldn't be fixable generally; can't you at least project the axis ranges just as you project the data values?

[philippjfr]

If you're referring to longitude being offset by 180 degrees, what's the underlying cause of that?

The actual projection is the cause of that, it's what cartopy outputs when projecting to crs.PlateCarree(central_longitude=180).

It seems strange that it wouldn't be fixable generally; can't you at least project the axis ranges just as you project the data values?

That is exactly what is being done, the axes simply reflect the data. The fix will be to make them not reflect the data, i.e. applying some custom formatting to offset the value.

[jbednar]

Ah, I see; this projection is centered with zero values near the international date line instead of Africa, so it's no longer correct to label the axis "Longitude". I wonder what people call that axis when using those numbers! Anyway, I agree, transforming the numbers back to Longitude would make the most sense to me, though I guess I still don't quite get why it couldn't just be done like that for all projections. No need to explain further, though; thanks!

[philippjfr]

though I guess I still don't quite get why it couldn't just be done like that for all projections. No need to explain further, though; thanks!

It's worth explaining, in the case of PlateCarree and RotatedPole the transform is simply translation of the xaxis, which we can easily implemented as a custom formatter. Other projections on the other hand require much more complex transforms to be mapped back to sensible latitude/longitude coordinates. This is still somewhat problematic because the lat/lon hover info on projections other than Mercator and PlateCarree will be wrong (the axes are already automatically hidden in these cases).

[jbednar]

Makes sense, thanks.