Raster colorization via "raster-color"

[rreusser]

This PR implements a raster-color style spec property as originally requested in #3889. It is closely analogous in both usage and implementation to heatmap-color. The essence of this feature is to reduce an RGBA pixel input to a scalar value, then to map that to a one-dimensional color scale. For efficiency, a lookup table is tabulated on the CPU and then sampled on the GPU, so that performance impact should be minimal.

This PR implements the following style-spec paint properties:

• raster-color (expression): defines a color map by which a raster layer is colorized. Must implement a style spec expression which operates on ["raster-value"] and returns a color.
• raster-color-mix (number[4]): specifies the combination of red, green, and blue channels used to reduce input to a scalar value for color map lookup. The scalar value is computed from the source color src using the formula: src.r * mix.r + src.g * mix.g + src.b * mix.b + mix.a. Default is luminosity, i.e. [ 0.2126, 0.7152, 0.0722, 0 ].
• raster-color-range (number[2]): defines the range in which the color map is tabulated. The raster-color expression is evaluated at 1024 steps (non-configurable) uniformly distributed in the range specified by raster-color-range.

The video below cycles through a number of use cases highlighted by the demo page, in particular:

1. Colorization of satellite imagery
2. Colorization of grayscale radar tiles (rreusser.b3c8wj9i, rreusser.5vzfcdjn)
3. Colorization of rgb-encoded terrain tiles (mapbox.terrain-rgb)
4. Colorization of rgb-encoded population density tiles (rreusser.796swou6)
5. Mapping of grayscale cloud layer (rreusser.dv2uv1dw) to transparent white layer (works with video layers as well, but omitted here for size; see, for example, the globe view cloud cover demo, which raster-colorization significantly simplifies)

Note: named color scales are not part of this feature. For a demo page, I've tabulated them using d3.

raster-color.mov

Hypsometric shading example

I think an example will make usage most clear. Let's consider what's basically the most nontrivial example, in which data is rgb-encoded. (RGB-encoding is not a requirement, and grayscale inputs permitting 256 levels, for example, require only specification of raster-color to style.) Terrain tiles are RGB-encoded via the formula height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1). To decode elevation, we can specify:

'raster-color-mix': [
  256 * 256 * 256 * 0.1,
  256 * 256 * 0.1,
  256 * 0.1,
  -10000
]

Note the extra factors of 256. This happens because [0, 255] is normalized to [0, 1] in the shader. We now have an elevation to work with, in meters. We can define the raster-color property:

'raster-color': [ "interpolate",
  ["linear"],
  ["raster-value"],
  0, "rgb(57, 143, 83)",
  500, "rgb(178, 205, 174)",
  1000, "rgb(221, 207, 153)",
  2000, "rgb(207, 155, 103)",
  4000, "rgb(227, 210, 197)",
  8848, "rgb(255, 255, 255)"
]

Note that the color map is parameterized in terms of meters, from sea level (0 meters) to mount everest (8848 meters). Then finally, we direct it to tabulate the color map in this range:

'raster-color-range': [0, 8848]

The final note is that RGB-encoded values must be sampled from uncompressed png tiles. For historical reasons, such tiles must be requested from Mapbox using the extension .pngraw. Furthermore, Mapbox may overwrite extensions when requesting tiles in order to minimize data transfer, so that you must use transformRequest to force requesting of PNG tiles. We must also avoid resampling, so that when requesting tiles from Mapbox we remove '@2x' as well. For example:

new mapboxgl.Map({
  ...,
  transformRequest: function (url, type) {
    if (type !== "Tile") return;
    return { url: url.replace('.webp', '.pngraw').replace('@2x', '') };
  }
})

The result looks something like this:

For many uses like colorization of grayscale imagery, the default values of raster-color-mix and raster-color-range would be fine and only raster-color would need to be specified.

Additional notes

• This comment pointed out that it's not pleasant to constantly transform between a normalized color map range (e.g [0, 1]) and physical values. This PR includes raster-color-range so that you can define the range over which the color map is defined and evaluated.
• Source alpha is not used in the raster-color-mix formula. Instead, the source alpha (that is, the alpha value of a tile image pixel) is carried through and applied as opacity to the final, colorized tile. The fourth component of raster-color-mix is instead always added as a constant offset. This feels a bit clumsy, but is what enables straightforward decoding of rgb-encoded data. For terrain, it would be the -10000 offset. See, for example, the above terrain example.
• Instead of first reducing to a scalar and colorizing in 1D, fully general 3D mapping of color, (r, g, b) -> (r', g', b'), was considered. However, implementation and usage are more complicated and the 3D color lookup table requires significantly more memory.
• This feature modifies the raster shader and so works wherever the raster shader is used. It therefore works with projections and video layers.

cc @mapbox/gl-native (It's my intent to work on implementation if/when we feel good about this PR. I believe all shader changes are contained by ifdef's so that this should be native-compatible, with the exception of me needing to figure out how to skip render tests.)
cc @mapbox/map-design-team for initial implementation

Launch Checklist

• briefly describe the changes in this PR
• include before/after visuals or gifs if this PR includes visual changes
• write tests for all new functionality
• document any changes to public APIs
• post benchmark scores
• manually test the debug page
• tagged @mapbox/map-design-team @mapbox/static-apis if this PR includes style spec API or visual changes
• tagged @mapbox/gl-native if this PR includes shader changes or needs a native port
• apply changelog label ('bug', 'feature', 'docs', etc) or use the label 'skip changelog'
• add an entry inside this element for inclusion in the mapbox-gl-js changelog: <changelog>Add raster colorization for raster layers</changelog>

[rreusser]

Included the sake of illustration and debugging, but we should probably remove tilesets on my account from the final PR.

[mourner]

This looks amazing, and such a cool feature for so little code! Especially love the API design and a big selection of demos. Great work @rreusser, let's get this in!

[endanke]

From gl-native perspective it looks good! There isn't too many changes and I can port it right away after it's merged.

[rreusser]

@endanke Thanks for the review! 🙇 Nothing is set in stone yet, so if you have any other opinions on the design, nothing is off the table. Points of contention/compromise, as I see it:

• Supports linear combinations only via raster-color-mix, so that in-shader log-scaling of data (or similar) is not supported
• The mapping only supports RGB -> scalar -> R'G'B'. @kkaefer suggested maybe a fully general RGB -> R'G'B' mapping using a 256x256x256 lookup table, but we opted for the simpler approach.
• The fourth component of raster-color-mix stands out in that it's just an offset which is not multiplied by source alpha. This makes decoding RGB-encoded data easy but feels just a bit low-level. I think it's okay though and just needs clear explanation or presets in studio. 🤷

Anyway, just wanted to get any decisions which I still find slightly unsettling on the table before etching it in stone. Thanks again!

[SnailBones]

Amazing feature @rreusser! 😍

One question: would it be possible to by default infer raster-color-range from the expression used in raster-color? Or would this add too much complexity/inconsistency with other expressions given the inner contents of one pre-evaluation expression would be influencing the behavior of another?

[rreusser]

@SnailBones I think it's possible, but the default value of raster-color-range would then have to be null instead of [0, 1] so that we could differentiate between the default [0, 1] as an affirmatively set value and as an unspecified value, which I guess would have to be applied later, where used. I get where you're coming from, but unless there's a good precedent for this, I think it would cause downstream uses like studio to have to know that there are two default values: on in the style spec, and one applied after the fact in software.

[rreusser]

One main thing has been holding this up: ensuring this will handle everything we want to throw at it. It's been suggested to use this for categorical data, so I've outlined below why this may be challenging.

How does categorical data differ from quantitative?

Categorical data means different input values are mapped to different colors, but that the order has no meaning. It's an encoding but not necessarily an ordering. Categories A, B, and C may correspond to RGB values 0, 127, and 255, but they could equally well be reordered. Quantitative data, of course can be mapped to 0-255 but the ordering of each value defines the meaning.

Why does this matter?

Colorization of categorical data results in artifacts at the boundary between categories.

The reason comes down to order of operations. There are two operations relevant here: 1) interpolation, and 2) color scale sampling. The order in which these are performed makes a big difference.

As a test, I made a simple texture with solid colors and no anti-aliasing and applied the d3 spectral color scale.

Desired output

We expect simple colorization, e.g.:

Actual output

What we actually get is boundaries that traverse the entire color scale between categories, rather than categories meeting smoothly.

Why does this happen? The average of black and white is 50% gray, so when we look up that color scale value, it could be anything. Below, input rgb(0,0,0) meets rgb(255,255,255). Due to linear interpolation, the result traverses the entire color scale and we get a whole rainbow where what we really want is direct interpolation from red to purple.

Where does this interpolation happen?

This interpolation happens in two places:

1. MTS. During tiling of GeoTIFF or similar input, colors are interpolated linearly. For this reason, it's avoided at a low level for historical terrain-rgb but can't be turned off for general tile sets.
2. Raster paint properties. Shader interpolation can simply be turned off with raster-resampling: 'raster-resampling': 'nearest'

So deactivation at the Mapbox GL layer level works fine if you have full control over the tiling. If you don't control tiling, then this problematic interpolation is baked into the tileset and there's no going back.

Can we work with this?

If we'd like the ability to display categorical data, then raster-sampling must be set to nearest and we must control the tiling process to avoid baking interpolation into the tileset. If those conditions are met, then it is possible to colorize categorical data.

[AlexanderBelokon]

Do we actually want direct interpolation from red to purple? Wouldn't it make more sense to have nearest neighbor or something like hqx for categorical data instead?

Also non-linear interpolations in MTS is a good feature to have in any case, since linear interpolation looks nasty, and the price for a nicer interpolation only has to be paid once. Is that on the roadmap at all? If it is, we could add nearest to the list of options too.

[rreusser]

@AlexanderBelokon That's a good point. The above presumes that nearest neighbor sampling is acceptable. Nicer sampling is, of course, very much preferable, but in order to accomplish that I think we'd have to (re)implement linear or better filtering in the shader, after colorization, as opposed to relying on built-in linear texture filtering.

[rreusser]

This PR will happen, but we're putting it temporarily on hold as we more carefully consider the API design before permanently committing to it.

[endanke]

@rreusser could we do some automatic switching between nearest and linear in the shader level? So we would have two color maps, one with nearest, one with linear. Then we would need to sample the neighbouring in the original image, and if the contrast between the neighbours is larger than a given threshold then we use the nearest variant, otherwise we use linear.
Not sure if it would totally solve the issues and of course selecting a threshold is also not obvious, but it might worth a try.

[pozdnyakov]

@rreusser should we provide a default value similarly to heatmap-color

      "default": [
        "interpolate",
        [
          "linear"
        ],
        [
          "heatmap-density"
        ],
        0,
        "rgba(0, 0, 255, 0)",
        0.1,
        "royalblue",
        0.3,
        "cyan",
        0.5,
        "lime",
        0.7,
        "yellow",
        1,
        "red"
      ],

[rreusser]

I don't have strong feelings on the default. I thought maybe grayscale was reasonable because we can really just say a lot less about the meaning of the data than we can for heatmap (which is necessarily point density).

[ghalter]

Hi, just wanted to ask if there is any news on this or an estimate on when this PR will be merged. We are really looking forward to this longed-for feature :-)

[mathewantony31]

@ghalter, thanks for your interest! We're still thinking through the details of this design and when it will become generally available. We're aiming to have more details on our plan by the summer (late June).

[shermify]

I love this feature and am probably using it a lot more than you intended at this point lol, but I just want to give some feedback from a future end-user. It works perfectly for applying color maps to grayscale imagery. Like you said there are some issues with categorical data, but generally resolved with 'nearest' resampling. However, that is a bit buggy on zooming in and out, etc. I would say categorical data is going to be a very common use case so probably worth getting right.

Super stoked for this feature, not sure how much I can contribute, but I have a bunch of complex use cases up and running if you need help testing. This is gonna be super popular when it hits.

[rreusser]

Thanks for the feedback, @shermify! It's very helpful! 🙇 In general I think there's a lot to be said for the "less is more" philosophy, that properties should be added sparingly and not every conceivable variant needs to be configurable. The categorical use case does keep arising though, and I think it mostly comes down to the distinction between colorize then interpolate (categorical) vs. interpolate then colorize (quantitative). I wonder if something like a "raster-color-mode": "quantitative" | "categorical" style spec property would be a sensible and necessary way to distinguish between the two.

[shermify]

I agree, simpler is better. I don't know if it even needs to be explicitly supported as a mode if it can work with tweaking existing settings. I have a map with lots of categorical data that needs to be mapped to many colors and you can see how the issue crops up.

However, setting raster-resampling to "nearest" works a bit better. You also need to set raster-fade-duration to 0 or it will "flash" the boundaries when tiles load in and it seems like there are some small rendering bugs when zooming in and out where the boundaries will suddenly pop in and "stick" for some reason. Other than that it seems to work fine. Performance is great even with a lot of categorical data and color.

[michaelbayday]

hey folks, us at cruise are greatly interested in this functionality and would appreciate a new version of mapbox to include this. i will also be filing up with a support ticket to understand when this will land

[shermify]

FYI for anyone following this issue, this feature is currently in Mapbox v3.0.0-beta. It was merged in via a different branch.

817fe63

[rreusser]

Thanks for pointing this out, @shermify! I'm indeed going to close this PR as the work made its way into the beta through a different commit. 🎉