Aggregating from non-orthogonal coordinate systems

[jonmmease]

Overview

This issue is to discuss the possibility of supporting 2D non-orthogonal coordinate systems in Datashader.

Current State

The Canvas constructor accepts x_axis_type and y_axis_type string arguments that may be set to 'linear' or 'log'. These strings are then used to lookup predefined axis objects, which are stored as the x_axis and y_axis properties of the Canvas instance:

https://github.com/pyviz/datashader/blob/3171d8869f63214cc161da755e88c30abe8b5098/datashader/core.py#L122-L147

The axis objects are subclasses of the datashader.core.Axis class:

https://github.com/pyviz/datashader/blob/3171d8869f63214cc161da755e88c30abe8b5098/datashader/core.py#L21

https://github.com/pyviz/datashader/blob/3171d8869f63214cc161da755e88c30abe8b5098/datashader/core.py#L90

https://github.com/pyviz/datashader/blob/3171d8869f63214cc161da755e88c30abe8b5098/datashader/core.py#L103

These classes have mapper and inverse_mapper methods that apply and reverse the 1D axis transform. e.g.:

https://github.com/pyviz/datashader/blob/3171d8869f63214cc161da755e88c30abe8b5098/datashader/core.py#L103-L114

These mapper and inverse_mapper methods must be numba jit compiled because they are called by the inner loops in the Glyph aggregation pathways. e.g.

https://github.com/pyviz/datashader/blob/3171d8869f63214cc161da755e88c30abe8b5098/datashader/glyphs.py#L158-L175

Current Limitations

1. It is not currently possible for a user to specify a custom axis type and pass that to the canvas constructor.
2. Because the Axis mapper and inverse_mapper interface only accepts a single value to map, these axes are fundamentally 1D. In other words, a transformed x value may depend only on the original x and it may not depend on original y value. This limitation precludes support for non-orthogonal 2D coordinate systems. e.g. polar, affine, barycentric, curvalinear, map projections etc.

Potential Enhancements

1. To allow users to customize the axis transformations, the datashader.core.Axis interface could be documented and the Canvas constructor could accept Axis sub-class instances as the x_axis_type and y_axis_type arguments.
2. To support non-orthogonal coordinate systems, I think it would be enough to update the mapper and inverse_mapper methods to accept two arguments.

def mapper(coord, other_coord):
    ...

def inverse_mapper(coord, other_coord):
    ...

Orthogonal transforms, like the current LinearAxis and LogAxis could ignore the second argument. Non-orthogonal transforms could make use of both arguments. Here's what a polar transform might look like.

class PolarRadiusToX(Axis):
    def mapper(r, theta):
        return r*np.cos(theta)

    def inverse_mapper(x, y):
        return np.sqrt(x**2 + y**2)

class PolarThetaToY(Axis)
    def mapper(theta, r):
        return r*np.sin(theta)

    def inverse_mapper(y, x):
        return np.atan2(y, x)

canvas = Canvas(..., x_axis_type=PolarRadiusToX(), y_axis_type=PolarThetaToY())
rs = ...
thetas = ...
canvas.points(x=rs, y=thetas)

Impact

This would require tracking down every use of Canvas.x_axis and Canvas.y_axis and adding the second argument, but I don't foresee it requiring any significant refactoring. This assumes that there won't be any performance degradation from adding the second argument to the Linear and Log axis mapping methods.

cc @jbednar @philippjfr

[jonmjoyce]

This is a feature I'm interested in. I am trying to rasterize weather model grids in an equirectangular projection but I was unable to find a way to do this with datashader. I see many examples with GeoViews but currently we just want to save images rather than integrate the whole holoviz stack.

[jbednar]

@jonmjoyce, I'd be happy to see what this issue proposes appear, but I don't think it's necessary for what you are asking for. It sounds like you already have data that can be accepted by the Datashader rectilinear quadmesh type, and if so just load that data into an xarray DataArray with independent x and y coordinates and call Canvas.quadmesh(da) to regrid it into a regular raster for displaying as an image. The result will effectively be a PlateCaree projection of the data, with each mesh cell scaled to cover the range of lat,lon values indicated by the equirectangular coordinates of that DataArray. If you want it to stay in the original coordinate system rather than mapping to PlateCarree you can replace those coordinates with some regularly spaced values (e.g. integers in order), in which case the output will be a straight scaleup or scaledown of the original data. And if you want some other (orthogonal) coordinate system, just project the coordinate arrays for x and y into some other projection (e.g. Web Mercator) before calling quadmesh.

We don't have any support for directly rendering into some irregular grid at the output side, but if you want is an image, I don't think you need that. And your data as described is already in an orthogonal coordinate system, which is also already supported (as indeed are curvilinear quadmeshes, just not curvilinear projections for points as this issue would address).

[jbednar]

See https://datashader.org/user_guide/Grids.html for examples of using grids like this in pure Datashader.