Datashader receives wrong range infromation from geoviews axis when viewing a cross section

[wx4stg]

versions

• bokeh 3.4.2
• cartopy 0.23.0
• datashader 0.16.3
• geoviews 1.12.0
• holoviews 1.19.1
• hvplot 0.10.0
• jupyterlab 4.2.3
• python 3.12.4
• Brave 1.67.123 Chromium: 126.0.6478.126 (Official Build) (64-bit)

I have a geoviews axis (lat and lon) of points that have multiple attributes attached to them, one of them is altitude. I want to have a cross section of longitude vs altitude above the plan axis. This works well.

However, when I enable datashader, the x-range of the plot does not update when I zoom in, so datashader does not return the correct pixel size of the rasterized images.

Here's a live example of what I'm talking about...
https://colab.research.google.com/drive/1sgYc2Y9HyEtMmggxMrRakzLw6Og53isv?usp=sharing

Code to reproduce

[lma.parquet](https://wx4stg.com/public/lma.parquet)

import pandas as pd

import holoviews as hv
import holoviews.operation.datashader
import datashader
import panel as pn

import geoviews as gv

hv.extension('bokeh')
gv.extension('bokeh')

ds = pd.read_parquet('lma.parquet')
easting3857, northing3857 = hv.util.transform.lon_lat_to_easting_northing(ds.lon, ds.lat)

data_plan = (easting3857, northing3857, ds.alt, ds.power)
kdims_plan = ['Longitude', 'Latitude']
vdims_plan = ['Altitude', 'Power']

data_cross = (easting3857, ds.alt, ds.lat, ds.power)
kdims_cross = ['Longitude', 'Altitude']
vdims_cross = ['Latitude', 'Power']

plan_points = hv.Points(data_plan, kdims=kdims_plan, vdims=vdims_plan)
plan_shaded = hv.operation.datashader.datashade(plan_points, aggregator=datashader.max('Power'), cmap='plasma').opts(
    width=300, height=300)
cross_points = hv.Points(data_cross, kdims=kdims_cross, vdims=vdims_cross)
cross_shaded = hv.operation.datashader.datashade(cross_points, aggregator=datashader.max('Power'), cmap='plasma').opts(
    width=300, height=100, ylim=(0, 20000))
pn.Column(cross_shaded, plan_shaded * gv.tile_sources.OSM())

Screenshots or screencasts of the bug in action

Screencast.from.2024-07-11.13-23-16.webm

[ahuang11]

Could be related to holoviz/geoviews#726

[ahuang11]

To reproduce easier:

import requests
from io import BytesIO

import pandas as pd

import holoviews as hv
import holoviews.operation.datashader
import datashader
import panel as pn

import geoviews as gv

hv.extension('bokeh', logo=False)
gv.extension('bokeh', logo=False)
data = BytesIO(requests.get('http://wx4stg.com/public/lma.parquet').content)
ds = pd.read_parquet(data)
easting3857, northing3857 = hv.util.transform.lon_lat_to_easting_northing(ds.lon, ds.lat)

data_plan = (easting3857, northing3857, ds.alt, ds.power)
kdims_plan = ['Longitude', 'Latitude']
vdims_plan = ['Altitude', 'Power']

data_cross = (easting3857, ds.alt, ds.lat, ds.power)
kdims_cross = ['Longitude', 'Altitude']
vdims_cross = ['Latitude', 'Power']

plan_points = hv.Points(data_plan, kdims=kdims_plan, vdims=vdims_plan)
plan_shaded = hv.operation.datashader.datashade(plan_points, aggregator=datashader.max('Power'), cmap='plasma').opts(
    width=300, height=300)
cross_points = hv.Points(data_cross, kdims=kdims_cross, vdims=vdims_cross)
cross_shaded = hv.operation.datashader.datashade(cross_points, aggregator=datashader.max('Power'), cmap='plasma').opts(
    width=300, height=100, ylim=(0, 20000))
pn.Column(cross_shaded, plan_shaded * gv.tile_sources.OSM())

[ahuang11]

I think to get around this, you can use geoviews for your map, and manually project for cross section; I use hvplot here because I was testing if I could get around it with hvplot, but it didn't work initially.

import requests
from io import BytesIO

import pandas as pd

import holoviews as hv
import geoviews as gv
import datashader
import hvplot.pandas
import cartopy.crs as ccrs
from holoviews.util.transform import lon_lat_to_easting_northing

hv.extension("bokeh", logo=False)
gv.extension("bokeh", logo=False)
data = BytesIO(requests.get("http://wx4stg.com/public/lma.parquet").content)
ds = pd.read_parquet(data)
x, y = lon_lat_to_easting_northing(ds["lon"], ds["lat"])
ds["x"] = x

(
    ds.hvplot(
        x="x",  # easting3857
        y="alt",
        kind="scatter",
        c="power",
        cmap="plasma",
        datashade=True,
        ylim=(0, 20000),
        frame_width=580,
        xaxis=False,
    ).redim.label(x="Longitude")   # match xdim of map
    + ds.hvplot(
        x="lon",
        y="lat",
        kind="points",
        c="power",
        cmap="plasma",
        datashade=True,
        tiles=True,
        crs=ccrs.PlateCarree(),
        projection=ccrs.GOOGLE_MERCATOR,
        frame_width=600,
    )
).cols(1)

[ahuang11]

Also, this issue probably is more holoviews/datashader related than geoviews since your example doesn't use geoviews at all.

Edit: actually it does, but you can call tiles from HoloViews:

hv.element.tiles.OSM()

[hoxbro]

Are you sure this is a HoloViews problem? Holoviews tiles are different from Geoviews tiles. Mainly, IIRC, no transformation is done with HoloViews tiles.

[ahuang11]

I mean it works with geoviews, but stripping out all the calls to geoviews, using vanilla HoloViews, the problem arises (demoed in the video)

import pandas as pd

import holoviews as hv
import holoviews.operation.datashader
import datashader
import panel as pn

hv.extension('bokeh')


data = BytesIO(requests.get("http://wx4stg.com/public/lma.parquet").content)
ds = pd.read_parquet(data)
easting3857, northing3857 = hv.util.transform.lon_lat_to_easting_northing(ds.lon, ds.lat)

data_plan = (easting3857, northing3857, ds.alt, ds.power)
kdims_plan = ['Longitude', 'Latitude']
vdims_plan = ['Altitude', 'Power']

data_cross = (easting3857, ds.alt, ds.lat, ds.power)
kdims_cross = ['Longitude', 'Altitude']
vdims_cross = ['Latitude', 'Power']

plan_points = hv.Points(data_plan, kdims=kdims_plan, vdims=vdims_plan)
plan_shaded = hv.operation.datashader.datashade(plan_points, aggregator=datashader.max('Power'), cmap='plasma').opts(
    width=300, height=300)
cross_points = hv.Points(data_cross, kdims=kdims_cross, vdims=vdims_cross)
cross_shaded = hv.operation.datashader.datashade(cross_points, aggregator=datashader.max('Power'), cmap='plasma').opts(
    width=300, height=100, ylim=(0, 20000))
pn.Column(cross_shaded, plan_shaded * hv.element.tiles.OSM())