Feature request: Bivariate colormapping

[sakunpanthi]

Wolfram recently launched Mathematica version 12 promising the support of complex plotting in 2d and 3d both.It is both wonderful and great for visualization! Matplotlib package also needs a color plotting system for complex functions!!

[story645]

Can you please provide more detail about what sort of features/plots you're proposing? If possible, a sample API would be extremely helpful, but even some examples of the types of plots you're asking for would be great.

[jklymak]

https://reference.wolfram.com/language/ref/ComplexPlot.html The big thing seems to be the ability to make the complex phase or amplitude have a shading. I actually think a shaded pcolor would be a useful addition in many realms.

shadedpcolor(x, y, Z1, Z2, norm1, norm2)

See also hill shading.. https://matplotlib.org/gallery/specialty_plots/topographic_hillshading.html

[QuLogic]

That sounds like #8738 (once it's complete)? Though I feel like I've seen an example of this sort of thing somewhere...

[u55]

I can already get something similar with alpha blending. Not quite perfect yet. I used seaborn for a nearly-constant-intensity cyclic colormap. To my eyes, the colormap seems to have some banding at the orange, green, and blue boundaries.

Compare with Mathematica.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, LogNorm
from matplotlib import cm
import seaborn as sns

# Cyclic colormap
my_cmap = ListedColormap(sns.color_palette("husl", 256))
# Alpha-blending grayscale
alpha_cmap = cm.gray(np.linspace(0, 1, 256))
alpha_cmap[:, -1] = 1 - np.kaiser(256, 15)
alpha_cmap = ListedColormap(alpha_cmap)

def f(x, y): # arbitrary complex function with 1 pole and 3 zeros
    z = x + 1j*y
    return z**2 - 8/z

x = y = np.linspace(-3, 3, 1000)
X, Y = np.meshgrid(x, y)
C = f(X, Y)

fig, ax = plt.subplots()
im = ax.imshow(np.angle(C)/np.pi, cmap=my_cmap, vmin=-1, vmax=1, extent=(-3, 3, -3, 3), origin='lower')
ax.imshow(np.abs(C), cmap=alpha_cmap, norm=LogNorm(0.1, 200), extent=(-3, 3, -3, 3), origin='lower')
cbar = plt.colorbar(im)
cbar.set_label('phase ($\pi$)')
ax.set_xlabel('real')
ax.set_ylabel('imaginary')
fig.tight_layout()
plt.show()

[tacaswell]

My first thought reading the issue title was "but the most common criticism we get is the plotting system is too complex!" 😜

[sakunpanthi]

Can you please provide more detail about what sort of features/plots you're proposing? If possible, a sample API would be extremely helpful, but even some examples of the types of plots you're asking for would be great.

Sorry I just couldn't attach it early! Here it is.

[sakunpanthi]

I can already get something similar with alpha blending. Not quite perfect yet. I used seaborn for a nearly-constant-intensity cyclic colormap. To my eyes, the colormap seems to have some banding at the orange, green, and blue boundaries.

Compare with Mathematica.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, LogNorm
from matplotlib import cm
import seaborn as sns

# Cyclic colormap
my_cmap = ListedColormap(sns.color_palette("husl", 256))
# Alpha-blending grayscale
alpha_cmap = cm.gray(np.linspace(0, 1, 256))
alpha_cmap[:, -1] = 1 - np.kaiser(256, 15)
alpha_cmap = ListedColormap(alpha_cmap)

def f(x, y): # arbitrary complex function with 1 pole and 3 zeros
    z = x + 1j*y
    return z**2 - 8/z

x = y = np.linspace(-3, 3, 1000)
X, Y = np.meshgrid(x, y)
C = f(X, Y)

fig, ax = plt.subplots()
im = ax.imshow(np.angle(C)/np.pi, cmap=my_cmap, vmin=-1, vmax=1, extent=(-3, 3, -3, 3), origin='lower')
ax.imshow(np.abs(C), cmap=alpha_cmap, norm=LogNorm(0.1, 200), extent=(-3, 3, -3, 3), origin='lower')
cbar = plt.colorbar(im)
cbar.set_label('phase ($\pi$)')
ax.set_xlabel('real')
ax.set_ylabel('imaginary')
fig.tight_layout()
plt.show()

Don't know if its going to work for complex functions with height and shading both as output.

[github-actions]

This issue has been marked "inactive" because it has been 365 days since the last comment. If this issue is still present in recent Matplotlib releases, or the feature request is still wanted, please leave a comment and this label will be removed. If there are no updates in another 30 days, this issue will be automatically closed, but you are free to re-open or create a new issue if needed. We value issue reports, and this procedure is meant to help us resurface and prioritize issues that have not been addressed yet, not make them disappear. Thanks for your help!

[jklymak]

I'll keep this (and changed the title) as I don't think bivariate colormapping is tracked in any open issues. #8738 was an attempt at this, but I deeply object to the proposed API there. but something like bivariate_pcolormesh and bivariate_colorbar would be useful and tractable projects.

[story645]

bivariate_pcolormesh

I'd be very against this b/c it would very quickly grow to bivariate_*{every method that takes a scalermappable} b/c theoretically every method that takes a univariate colormap should be able to scale up to nvariate. This also sorta ties into the alpha_cmap discussion - changing alpha based on the data...and I think the solution may be more around a more customizable/composable approach to colormaps.

[jklymak]

Which other methods do you think would look good with a bivariate colormap?

[story645]

Which other methods do you think would look good with a bivariate colormap?

Honestly look good is besides the point - I think the library should facilitate creating visualizations in a consistent matter. Bi/nvariate colormaps are valid colormaps and therefore should be allowed everywhere where we allow univariate colormaps. Here's an implementation that seems to do that xycmap

If I had to start somewhere, one of the most common use cases is for choropleths, which means patchcollection needs to make colorbars. Here's a pysal example.

Also scatter plots/multivariate plots because there's a sort of gotta encode em all in viz where folks will try to use whatever tools are - available to encode variables. Also heatmap for the discrete discontinuous case.

If I'm trying to make a bivariate heatmap, I'm probably not distinguishing between imshow and pcolor unless I have to.

Which is back to my original point - if at all possible, I'd prefer the library to be as unopinionated as possible on where folks can use specific types of colormaps b/c we don't limit any of the other ones.

[jklymak]

A bivariate pcolormesh (or imshow, or chlorpleths) does not simply use a specific type of colormap - it requires two inputs, one for each of the variables, and it requires two colormaps. Thats very different from single mapping, so I don't see how you will cleanly snake the extra information through the existing API of all the methods.

The API choices are 1) try to guess that the input data are shaped in such a way that the user wants a bivariate colormap, 2) add extra kwargs to all the methods or 3) just have new methods.

I personally feel 3) is the more practical route; most of our colormappables take RGBA, so any new methods we write would just be wrappers that create the RGBA matrices to pass to the underlying methods. Having bivariate_pcolormesh is a straight forward and practical way to get the feature, compared to somehow snaking pcolormesh([A, B], cmap=['Reds', 'Greys'], clim=([limA, limB])) through the existing API. That will be a messy process, and all for a relatively obscure feature.

As for colorbars, there are probably going to need to be a few options. For a shaded colormapping, the colorbar can almost be the same as the default - the shading is obvious. For a more complicated blend, you'd want something like what xycmap is doing. Again, I think this would be better as separate methods because it requires different knobs to get the desired effect.

[story645]

1. try to guess that the input data are shaped in such a way that the user wants a bivariate colormap

I don't see why we can't make this explicit? If they are using a bivariate colormap then they must pass in a bivariate norm and then the must pass in data in a shape that the norm can take as input to make a color. Yes it changes the pipeline but it would in general make colormapping much more flexible if we tied the requirements of the input being mapped to color to the function that changes it into color.

Generalizing the colormapping pipeline or creating a new family of functions are both gonna have maintenance burdens; I think the former will require more investment upfront but that long term would be more consistent with the current library architecture that treats data->color as somewhat independent of the plot type.

I agree that for colorbars, we'd probably be better off with a .colorbar_nvariate(dims=n) then trying to shoehorn all that into the existing colorbar.

[jklymak]

I don't see why we can't make this explicit? If they are using a bivariate colormap then they must pass in a bivariate norm and then the must pass in data in a shape that the norm can take as input to make a color.

To me, that sounds worse for the user:

cmap = bivariate_cmap(cmap1, cmap2)
norm = bivariate_norm(norm1, norm2)
data = np.stack([data1, data1])
ax.pcolormesh(data, norm=norm, cmap=cmap, kwargs=kwargs)

plus all the maintenance burden of snaking that all through pcolormesh.

Conversely:

ax.bivariate_pcolormesh(data1, data2, norms=(norm1, norm2), cmaps=(cmap1, cmap2), kwargs=kwargs)

which then just translates to an RGBA array that gets passed to the usual pcolormesh.

[story645]

To me, that sounds worse for the user:

cmap = bivariate_cmap(cmap1, cmap2)
norm = bivariate_norm(norm1, norm2)
ndata = np.stack([data1, data1])
ax.pcolormesh(data, norm=norm, cmap=cmap, kwargs=kwargs)

This seems great to me: it's consistent with the ways we already make users work with specialized norms and colormaps & it's far easier to teach and publicize "hey we now support bivariate cmaps and here's how you use it" then to have to roll out a new bivariate_* function every release cycle or two.

I think your argument makes a ton of sense if the only function w/bivariate support we'd want to implement is pcolormesh. But I think having all the functions that support colormapping support nvariate colormapping is a reasonable request, and at that point I think it makes more sense to try to build a consistent pipeline we can try to reuse across functions than writing a new set that will probably each end up doing their own things in ways that will lead to inconsistencies - units are the prime example of where this is already happening.

[jklymak]

units are the prime example of where this is already happening.

Units are the prime example of why I think this is a bad idea! Our unit support is a brittle, inconsistent mess.

[story645]

Units are the prime example of why I think this is a bad idea! Our unit support is a brittle, inconsistent mess.

Yes because we do it on a per function basis, which is I think what you're proposing here for the special case of bivariate norm. I'm saying if you want consistency, then there should be a shared colormapping/norming pipeline.

[jklymak]

Yes because we do it on a per function basis, which is I think what you're proposing here for the special case of bivariate norm.

I don't think that is the case. We have a general units pipeline that is intricate, and gets into trouble when used in the individual methods. You are proposing basically the same thing - somehow matching up data shape with norm class and then having each method deal with that.

At least units, however, have a claim to widespread use - bivariate color mapping is pretty specialized and I don't think is worth either the complexity or the effort of creating a whole integrated pipeline that sits inside the code base. xycmap demonstrates that this doesn't even have to be part of the core library to work fine.

[tacaswell]

https://github.com/matplotlib/data-prototype/blob/main/examples/mulivariate_cmap.py is a functioning example of this in the data-prototype. Handling this sort of thing and handling units uniformly boil down to solving the same problem.

it requires two inputs, one for each of the variables, and it requires two colormaps

I do not think this is correct, in the example above it is using two degrees of freedom in the input data (e.g. it is a 2D array of R^2) to map to RGB via HSV which mixes the input values per pixel in a non-separable way. Similarly, for the norm the options are to map to the unit disk (sphere) or square (cube), the first of which is not separable into two norms. I have used this scheme for actual science to show a very dense orientation field.

I agree with @story645 that we do not want to pick up a second version of everything that creates subclasses of ScalarMappale and instead want to generalize the API we already have.

I think we are already doing option 1 (guess based on the shape of the input). I'm going to focus on imshow but all of this applies to anything where we have one artist that internally maps from R^n user data to R^4 color (per pixel, per polygon, per marker, per whatever).

If we look at the norm -> colormap chain, that in R^1 -> R^4 (I really need to go back to #18487) in the case where the user passes an array with shape (N, M). However, we can also view this is the user passing use (N, M, 1) but we are forgiving on then input. Thus, ScalarMapple handles the uni-variate case. We also let the user pass us (N, M, 3) and (N, M. 4) arrays which we, at an implementation level, just pass through to the backends (which require an RGBA array (or just an RGBA color)). However, we can also understand this as supporting the 3-variate and 4-variate cases in only the special case where the norm is to the unit cube/hypercube and both the norm and the colormap are effectively the identity. We can then understand the current implementation as a performance driven fast-path.

I suspect it is the case that in a lot of the code we did bake in the R^1 a bit too hard (e.g. the resampling code in imshow), but in other cases things will "just work" as (if we can escape the special casing)

matplotlib/lib/matplotlib/cm.py

Lines 505 to 509 in f588d2b

	# This is the normal case, mapping a scalar array:
	x = ma.asarray(x)
	if norm:
	x = self.norm(x)
	rgba = self.cmap(x, alpha=alpha, bytes=bytes)

is naive to the dimensionality of both the input and intermediate data!

---

There are obviously a lot of API work to be done (e.g. norm.vmin and norm.vmax make no sense in higher dimensions, cube vs sphere like norms (no difference in 1D), what do we do about over/under in higher dimensions (do we only have inbounds / out of bounds), can we do run-time checks before draw-time to make sure the pipeline will work, ...), but it is pretty clear to me the current set of norms / colormaps and the R^1 or RGB(A) are special cases of something more general.

[jklymak]

@tacaswell that sounds fine in theory. But your R^2 data is just as easily written as a list or tuple of two R^1 objects. I'll argue that 99.99% of the time people will not be storing data as R^2, but instead as two R^1 objects (eg ds.speed, ds.direction, or amplitude and phase from the OP). So now to use this hypothetical pipeline we are asking users to unnaturally create an R^2 object.

Further, we have already blocked out R^3 and R^4 for RGB(A). So, are we going to skip those, but allow R^1, R^2, R^5, R^6...? That doesn't seem worth the abstraction since R^2 is hard enough to make sense of, I doubt there is much call for higher order dimensions. Or are we going to somehow have an IdentityNorm that allows the R^3 and R^4 to pass through?

There are obviously a lot of API work to be done

Agreed ;-). Looking at all the changes in #8738 it is certainly non-trivial.

Look at it this way - adding bivariate_pcolormesh (or bivariate_*) doesn't at all preclude the pipeline you are proposing. If the new pipeline is sorted out, bivariate_*-type methods just become vestiges of a simpler time, like plot_date currently is.

BTW just to be clear, I too have used bivariate color mapping for real science: https://tos.org/oceanography/assets/docs/25-2_klymak.pdf

[tacaswell]

I'll argue that 99.99% of the time people will not be storing data as R^2, but instead as two R^1 objects

That has not been my experience; I have mostly dealt with vectors in Cartesian rather than radial representation (but I mostly look at small (effectively) flat things not the surface of a very big round thing ;) ).

In the general case where the norm is some sort of coupled ufunc we are going to have to make stacked numpy array eventually...but we are pretty good about being able to do that, imshow([dir, phase],...) should presumably work.

So, are we going to skip those, but allow R^1, R^2, R^5, R^6...?

I was thinking the default norm/colormap at R^3 and R^4 would trigger the pass-through logic (if isinsance(norm, IdentiyNorm): is not any worse than if data.ndim == 3:).

An overlay of elemental composition setting the color (element per RGB), alpha set by the sum of those elements on top of a gray-scale image of the sample via another mechanism is something that some of the beamlines at NSLS-II use and I think would be made easier by this sort of scheme.

Look at it this way - adding bivariate_pcolormesh (or bivariate_) doesn't at all preclude the pipeline you are proposing. If the new pipeline is sorted out, bivariate_-type methods just become vestiges of a simpler time, like plot_date currently is.

However I think that if we sort out the API issues for the R^2 case we will have solved them from the R^N case as well.

That said, I see how to do R^2 legends...I do not see how to do higher !

[trygvrad]

I will, but it will probably take some weeks, as I have to prioritize my day job for a while.
I'll try to keep you updated

[story645]

Awesome! And please keep in mind that you can always open up a draft PR to show sketches of things

[trygvrad]

Short update:
I'm prototyping this feature in https://github.com/trygvrad/matplotlib/tree/mutivariate_colormaps
And building colormaps in https://github.com/trygvrad/multivariate_colormaps/
The current version allows:

...
fig, axes = plt.subplots(1, 2, figsize = (12,4))
pm = axes[0].pcolormesh((A,B,C))
fig.colorbar(pm)
pm = axes[1].pcolormesh((A,B,C), cmap = '3VarSubC', vmax = (0.4, 0.7, 0.5))
c0, c1, c2 = fig.colorbar(pm)

And all existing tests pass. (data from https://arxiv.org/abs/1812.10366)

I think this is the syntax we want, but there are details in the implementation that are worth discussing.
I will try to open a draft PR soon, so we can discuss the way cm.ScalarMappable has been adapted to accommodate multivariate data.
Even when we can agree on the syntax and architecture, a lot of work remains with implementing this for all the other visualization forms (I prototyped with pcolormesh) building dedicated tests, updating documentation etc.

[story645]

Awesome! We have new contributor calls first Tuesday of every month and weekly project calls every thursday. Please join if you can make either: https://scientific-python.org/calendars/

[jklymak]

I'll reiterate that I am deeply skeptical about overloading like this.

A) colorbar is pretty finicky as it is. Do all the kwargs work with this? What does shrink=0.6 do? Location=bottom? How does all the anchoring work? What happens when constrained layout is called?

B) I'd suggest drafting the new docstrings. Suspect this is going to be very hard to explain in a way that doesn't detract from the existing docs.

C) error checking is going to be a chore. What happens if I only specify two vmax values? What about only one?

So while I'm sure prototyping this is straightforward it is going to be a lot of code interwoven with our already complex code. We need to decide if this is worth the cost.

[trygvrad]

I'll try to make the new contributor call :)

A) I agree that colorbar is somewhat finicky. I built the current implementation by first designing a test that checks the various placements, shrink, etc., although it will need to be updated later to check constrained layout, see https://github.com/trygvrad/matplotlib/blob/mutivariate_colormaps/lib/matplotlib/tests/baseline_images/test_multivariate/multivariate_colorbars.png .
I find it intuitive to let the new colorbars be placed in the space that would otherwise be used by a single colorbar. You are free to propose a different behavior, but I ask that you remain open to us working together to find a solution.

B) I'd suggest drafting the new docstrings I already am :) But I believe that a collaborative approach greatly improves the final outcome. For pcolormesh the docsting currently includes the following, but I would be happy to hear your suggestions as to how this can be improved, and if/how the rest of the docstring should be adapted.

        Parameters
        ----------
        C : array-like
            The mesh data. Supported array shapes are:

            - (M, N) or M*N: a mesh with scalar data. The values are mapped to
              colors using normalization and a colormap. See parameters *norm*,
              *cmap*, *vmin*, *vmax*.
            - (M, N, 3): an image with RGB values (0-1 float or 0-255 int).
            - (M, N, 4): an image with RGBA values (0-1 float or 0-255 int),
              i.e. including transparency.
            - (V, M, N) for N > 4: V images with multivariate scalar data, in 
              this case *cmap*, *vmin*, *vmax* accepts lists of length V.

            (M, N) define the rows and columns of the mesh data.

C) error checking is going to be a chore. What happens if I only specify two vmax values? What about only one? I do not think these questions are as difficult as you make them sound. I already built a function that sanitizes the input and raises an error if the user input is ambiguous :

    def sanitize_multivariate_data(self, data, norm, cmap,
                                        vmin, vmax, multivar_mode):
        """
        If data is not multivariate, this function returns the input unchanged.
        If data is multivariate, i.e. it has shape (n, x, y) with 1<n<9 this
        function will ensure that norm, cmap, vmin and vmax are also of length
        n.
        
        - A single argument for norm, vmin or vmax will be repeated n times in the
        output.
        - If cmap is a string corresponding to a family of multivariate colormaps,
        the family will be returned as a list.
        - If multivar_mode is None, this is set to 'Add' unless the name of the first
        colorbar includes the substring 'Sub'.
        
        returns:
            n, norms, cmaps, vmins, vmaxs, multivar_mode
        """

The wording could probably do with some collaborative improvements, but the logic and error handling is there.

[jklymak]

I'd recommend the proper Thursday developers call for this.

    - (M, N, 3): an image with RGB values (0-1 float or 0-255 int).

        - (M, N, 4): an image with RGBA values (0-1 float or 0-255 int),
          i.e. including transparency.
        - (V, M, N) for N > 4: V images with multivariate scalar data, in 
          this case *cmap*, *vmin*, *vmax* accepts lists of length V.

Right, so this is already ambiguous if a user passes (3,3,3) what they meant.

I think you should try and get feedback from other senior devs: @timhoffm @QuLogic @efiring @anntzer @greglucas come to mind, apologies fir anyone else has worked on colorbar things that I've forgotten.

[jklymak]

Also @ksunden if for a typing point of view if nothing else.

[anntzer]

Not that I looked in depth at this at all, but (V, M, N) seems a bit unusual; more standard ordering would be (M, N, V)? (Notwithstanding the ambiguities that can occur in both cases.)

[trygvrad]

Regarding (N,M,V) vs (V,N,M) I think the 'default' data structure will depend largely on the data source.
That said, I would prefer (V,N,M) because it allows for the same input structure for data as for vmin, vmax, norm, etc. (at least in the case where the different variates are represented by different variables):
i.e.
(V,N,M) → pcolormesh((A,B,C), vmax = (0.4, 0.7, 0.5))
vs
(N,M,V) → pcolormesh(np.stack((A,B,C), axis = -1), vmax = (0.4, 0.7, 0.5))

[trygvrad]

I think we should also consider an alternative option: To differentiate between 3D arrays and an iterable of 2D arrays.
i.e.
np.empty(3,3,3) → RGB image
[np.empty(3,3), np.empty(3,3), np.empty(3,3)] → multivariate

Multivariate data will always be divided into separate 2D arrays for normalization, so there is no additional overhead on the backend if starting from a list.

[anntzer]

Diverging behavior between list-of-nd-arrays and (n+1)d-arrays already exists e.g. in boxplot() and I would personally consider that API to be a nightmare and not something to be expanded.

[jklymak]

I'll point out the other obvious issue here, which is that this API does not admit bivariate color mapping, but rather additive color mapping.

[trygvrad]

As I argued previously, I think it is sensible to build multivariate colormaps first, and then implement 2D colormaps.
Unless you wish to argue that we should not support multivariate colormaps at all?

May I ask: when you think about bivariate colormaps do you think about something like the figure on the left, or something like the figure on the right?

On an implementation side I think the two sides (left/right) in the figure above are quite different. It would seem sensible to me to first implement then functionality on the left (multiple 1D colormaps) in a generalized way, and then move on to develop the functionality on the right (a single 2D colormap). i.e. these are two different features, and I think they belong in separate pull requests (where the second one will partially build on the first).

[jklymak]

Yes, but the proposed API can only be used for one of the two.

Unless you wish to argue that we should not support multivariate colormaps at all?

My argument is that we should support them, but that they should have stand-alone methods, not overload the existing methods.

[jklymak]

Discussed on the call today. Next-steps were for @trygvrad to prototype VectorMappable (name up for discussion) and maybe write up or prototype how that would be threaded through the library. There are issues to do with documentation, and there are issues to do with differentiating between the "additive" color mapping and fully-nonlinear multivariate color mapping.

Colorbars were also briefly discussed and the leading idea is that these probably need standalone methods.

[trygvrad]

I won't make it to the meeting today, but I am currently building a VectorMappable in a new fork.
You can look at the diff to get an overview of the required changes.

In this version the cmap acts as a switch, if a string, the cmap is converted to BivarColormap or MultivarColormap if the string is in the corresponding database.
Type checks are then used to select the appropriate behavior.
I hope to make a draft PR before the meeting next week, and then we can discuss if this is a good approach at that time.

[trygvrad]

I have submitted a draft pr #26996
Ideally we can discuss this at the meeting tomorrow, and next week

[endolith]

#4369 Wishlist: Multidimensional colormaps

[scottshambaugh]

Wanted to highlight another implementation of 2D colormapping that has been bolted onto matplotlib by @fzaverl here: https://github.com/fzaverl/s3dlib#2d-colormaps

They did a really nice use case for 3D surfaces, where the u, v surface normals are taken as the two dimensions of the color map. More generally known as a normal map:

[trygvrad]

I still intend to see this through, but have been otherwise occupied for a while.

Besides the technical implementation, we need to find and/or design some bivariate colormaps.
I have been thinking about this for a long time, and designed a number of bivariate colormaps when I built my previous package colorstamps.
I finally got around to writing down the design principles I have been contemplating in a blog post.

If there is any other work of this nature that we should consider, it would be good to start collecting it now.

[endolith]

If there is any other work of this nature that we should consider, it would be good to start collecting it now.

Here's mine: https://github.com/endolith/complex_colormap

[story645]

we need to find and/or design some bivariate colormaps.

Preference is for colormaps backed by either an academic paper or reputable white paper (while also still being public domain)

[trygvrad]

we need to find and/or design some bivariate colormaps.

Preference is for colormaps backed by either an academic paper or reputable white paper (while also still being public domain)

I completely agree.
I find that this paper : Bernard, Jürgen, et al. “A survey and task-based quality assessment of static 2D colormaps.” Visualization and Data Analysis 2015. Vol. 9397. SPIE, 2015.
Provides a good overview of published 2D colormaps at that time (2015), and how they compare.
I believe some progress has been made in colormap design since then, with packages such as colorspacious (V1.0.0 Nov 2015) making design and evaluation of colormaps much easier.
When I look for more recent work, I find many that make use of these capabilities (such as @endolith), but I have yet to find any that discuss the broader context, such as how diverging/scalar/cyclic data needs different colormaps, or how multivariate colormaps intersects with colorblindness.

I would like to write a paper on this myself (ideally as a collaboration with other interested parties), but that may take a while, which is why I started with a blog post :)

On the other hand, I would be very happy if any of you know some more recent works that I have missed.