Using datashade=True flips the image up-down

[bidhya]

Description of expected behavior and the observed behavior

Using datashade=True option the plotted image is flipped up-down. Is there some bug due to datashade or I am doing something wrong. I also see the dynamic range of color changed as well.

See the two images below side-by-side, with only change being 'datashade=True'

Here dem is any 1-band DEM

Complete, minimal, self-contained example code that reproduces the issue

dem = dems[0] 
ds =rioxarray.open_rasterio(f'{base_folder}/{dem}')
print(f'DEM: {dem}')
print(ds.shape)

asp = ds.shape[-1]/ds.shape[1]

ds.hvplot.image(x='x', y='y', aspect=asp, frame_width=250, cmap='gist_earth') + ds.hvplot.image(datashade=True, x='x', y='y', aspect=asp, frame_width=250, cmap='gist_earth')

Screenshots or screencasts of the bug in action

ALL software version info

all these libraies from conda-forge
hvplot: 0.7.0
datashader 0.12.0
holoviews 1.14.1
xarray 0.16.2
rioxarray 0.1.1

[jbednar]

That's clearly a bug, thanks! There's an inverting-axis bugfix in panel 0.10.3 that was just released; can you update with conda install -c pyviz panel=0.10.3 and try your example again? I can't tell from the code above whether this specific fix applies here.

[bidhya]

I updated to panel 0.10.3 from conda-forge. Still the same issue.

Just out of curiosity, this bug should be in hvplot, holoviews, or datashader so why would updating panel package fix this issue.

[jbednar]

hvPlot is based on HoloViews, and HoloViews uses Panel to lay out figures, including making all the axes match and be normalized together. And Panel did have a bug where such code wasn't working with inverted axes, but if that fix doesn't help here, then it's probably an issue in HoloViews (though it could also be Datashader). Can you provide a reproducible test case?

[jlstevens]

This seems fairly important to me so I'll assign this to the milestone.

[bidhya]

Here is a minimum code to reproduce the issue. Github does not allow me to upload tiff image, so I will just post the code. But replace the example.tif with any geotiff

`
import xarray as xr
import hvplot.xarray
ds = xr.open_rasterio('example.tif')

ds.hvplot.image(rasterize=True) # 1. Using rasterize option
ds.hvplot.image(datashade=True) # Using datashade option, the axis is flipped
`

[michaelaye]

gosh, this is a biggie. It's actually not simply flipping the axis and the data, it flips only the shown data without flipping the axis, so the data alignment is lost and wrong axis values are shown for the pixels.
This is highly confusing when exploring data for which I have no expectation on how things look.

This has another annoying side effect that when zooming in, the delta-y used internally for the zoom operation placement does things correctly, only to be flipped by the datashade operation, IOW, the end result jumps wildly around in the opposite y-direction. This can nicely be seen by carefully do a pan in one y direction and seeing how the resulting data moves into the opposite direction. (Or, in my below linked MCVE notebook, one has to zoom into the upper left corner to end up showing the anomaly, despite being shown in the lower left. Zooming in on the anomaly gets you nowhere in the datashaded plot).

The issue only appears when the coordinate in question has descending values, which is very often the case for geospatial analysis, which is why I'm kinda puzzled why nobody complained before? Must be hitting hundreds of users at least unless they never use datashade?

I have designed an MCVE with comments here:
https://gist.github.com/3b1b4b6e713c49d6d1d021d8cbd13592

and for your convenience the minimal code until exposure of the basic issue pasted here:

import numpy as np
import xarray as xr
import hvplot.xarray
from holoviews import opts

arr = np.arange(20000).reshape(200, 100)
da = xr.DataArray(
    data=arr,
    dims=['x','y'],
    coords = {
        'x':np.arange(200),
        'y':np.arange(100,0, -1)
    }
)
# anomaly placed in UPPER left corner (due to descending y-coord)
da[10:20, 10:20] = 0
# anomaly shows up at wrong coordinates in LOWER left corner
da.hvplot(datashade=True, x='x', y='y', data_aspect=1)

[michaelaye]

Workaround:

Flip the descending axes with np.flip(obj, axis).
It will flip the data AND the appropriate coordinate accordingly so that datashader has no problem.
The default plotter plots ascending y in any case, so one can flip plots in the plot command if required.

For above MVCE:

da = np.flip(da, 1)  # for axis 1 where y is located

[jbednar]

We've definitely run into flipping issues in various cases before, but not this particular one; the files we test with must all have ascending coordinates. It would be great to have a check to catch this condition and do the flipping if needed.

[hoxbro]

By using @michaelaye MCVE as a starting point and some heavy debugging. I found the bug happens when applying holoviews.operation.datashader.shade to a holoviews.Image.

import holoviews as hv
import numpy as np
import xarray as xr
from holoviews.operation.datashader import shade
hv.extension("bokeh")

arr = np.arange(20000, dtype=np.float64).reshape(200, 100)
da = xr.DataArray(
    data=arr, dims=["x", "y"], coords={"x": np.arange(200), "y": np.arange(100, 0, -1)}
)
da[10:20, 10:20] = 0
image = hv.Image(da)
image + shade(image)

With this information it was pretty straight forward to find that the problem happens when calling datashader.utils.orient_array in

datashader/datashader/transfer_functions/__init__.py

Line 212 in 71ba1d5

data = orient_array(agg)

If I replace the function call with data = agg.data I get the same output with and without the shade.

import panel as pn
pn.extension()


def create_plot(x_reverse=False, y_reverse=False):
    x, y = np.arange(200), np.arange(100)
    if x_reverse:
        x = x[::-1]
    if y_reverse:
        y = y[::-1]

    da = xr.DataArray(
        data=arr, dims=["x", "y"], coords={"x": x, "y": y}
    )
    da[10:20, 10:20] = 0
    image = hv.Image(da)

    return pn.Column(f"# x reverse = {x_reverse}", f"#y reverse = {y_reverse}", image, shade(image))


pn.Row(create_plot(), create_plot(True), create_plot(False, True), create_plot(True, True))

Before:

After:

I could not determine if the fix is to "only" remove the function call to orient_array or if the orient_array should give the correct output. Nonetheless, I feel like this is a very good start to fixing the bug!

Note if you want the correct direction right now you can monkey patch the orient_array like this

import datashader
datashader.transfer_functions.orient_array = lambda x: x.data

[jlstevens]

@philippjfr

@MridulS and I agree this is one of the more severe bugs we may want to try to fix for 0.8? Do you agree?

[philippjfr]

As @hoxbro's investigation showed I believe this is extremely unlikely to be an hvPlot or even HoloViews bug. I do believe this is down to datashader.util.orient_array and/or how it's used by the shade function and by the Canvas.raster method in datashader.

[ianthomas23]

This is indeed a bug in datashader. It is simple to understand, but is going to be awkward to fix without breaking other use cases.

The clue is that GeoTIFF files usually (always?) have their dimensions in the order (y, x) and the y-coordinates are decreasing. Other uses still use the numpy standard of (y, x) but coordinates increasing, plus it is entirely possible to use (x, y) dimensions instead. In datashader/utils.py we have inconsistent handling of the dimension order:

• calc_res assumes (y, x) and flips the y-axis.
• calc_bbox assumes (x, y).
• orient_array assumes (y, x) as most of the work is done by calc_res.
• downsample_aggregate uses (x, y).

We could do a localised fix in calc_res to check the dimension order and behave correctly, but then I suspect that (possibly incorrect) assumptions made elsewhere in datashader and possibly even in HoloViews/hvPlot might be affected.

Various functions calculate and/or pass around array properties named similar to (xprop, yprop). In some cases these might actually be referring to properties in the (x, y) order and in others they might be referring to properties in the same order as the array indices.

[jbednar]

Fixed in #1095.