How to plot an observer on a georeferenced tif file

Datashader
1

Hello,
I am new to the forum, but I thought here would be the best place to ask my question, since I could not find the anser so far anywhere else.

I want to visualize a georeferenced tif and plot an observer as a dot on the canvas.

Here is a link to the tif file: https://www.dropbox.com/s/kmh68glncnla06u/dsm_sampled.tif

I have written the following code, based on a notebook I found on geoviews:

import datashader as ds 
import geoviews as gv
import xarray as xr
import pandas as pd
import datashader.geo as dsgeo
from datashader.transfer_functions import shade, stack,dynspread
from datashader.colors import Elevation

# read tif
tif = xr.open_rasterio("dsm_sampled.tif")
# take a band and replace nan values
terrain = tif[0].fillna(0)

# create canvas
width = 1000
height = 1000

cvs = ds.Canvas(plot_width=width, plot_height=height)
illuminated = dsgeo.hillshade(terrain)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')


# some arbitrary position
x_pos = -243793.842
y_pos = -252500

observer_df = pd.DataFrame({'x': [x_pos],'y': [y_pos]})
observer_agg = cvs.points(observer_df, 'x', 'y')

observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),
                            threshold=1, max_px=4)

####
# it is working if I leave out observer_shaded 

# stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),
#       terrain_shaded)
###

# This call is not working for some reason
stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),
      terrain_shaded, observer_shaded)`

But I am getting the following error and can not find a solution:

    ---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-ebf9dfe88990> in <module>
     40 # This call is not working for some reason
     41 stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),
---> 42       terrain_shaded, observer_shaded)

~/opt/anaconda3/envs/grass/lib/python3.7/site-packages/datashader/transfer_functions/__init__.py in stack(*imgs, **kwargs)
    129     imgs = xr.align(*imgs, copy=False, join='outer')
    130     with np.errstate(divide='ignore', invalid='ignore'):
--> 131         out = tz.reduce(tz.flip(op), [i.data for i in imgs])
    132     return Image(out, coords=imgs[0].coords, dims=imgs[0].dims, name=name)
    133 

~/opt/anaconda3/envs/grass/lib/python3.7/site-packages/toolz/functoolz.py in __call__(self, *args, **kwargs)
    301     def __call__(self, *args, **kwargs):
    302         try:
--> 303             return self._partial(*args, **kwargs)
    304         except TypeError as exc:
    305             if self._should_curry(args, kwargs, exc):

~/opt/anaconda3/envs/grass/lib/python3.7/site-packages/toolz/functoolz.py in flip(func, a, b)
    735     [1, 2, 3]
    736     """
--> 737     return func(b, a)
    738 
    739 

TypeError: ufunc 'over' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

If I am leaving out the observer point from the stack, everything is plotting fine. But with the observer point the plot breaks.

2

Sorry for the very late reply here. The problem (which should be detected with a better error message) is that your arrays do not line up with each other, i.e. without specifying an x_range and y_range on the Canvas it will fall back to auto-ranging and your single point will have a very different range than the other cmponents. You can fix this by computing the x_range and y_range from your terrain tif:

x_range = terrain.x.min().item(), terrain.x.max().item()
y_range = terrain.y.min().item(), terrain.y.max().item()

cvs = ds.Canvas(plot_width=width, plot_height=height, x_range=x_range, y_range=y_range)