How to use geoserver SLD style to serve single channel elevation raster ("gray" channel) as Mapbox Terrain-RGB tiles

Question

I have an elevation raster layer in my GeoServer with a single channel ("gray"). The "gray" values is elevations values (signed int16). I have 2 clients:

1. The first one is using that elevation values as is.
2. The second one expect to get [Mapbox Terrain-RGB format][1]

I do not want to convert the "gray scale" format to Mapbox Terrain-RGB format and hold duplicate data in the GeoServer. I was thinking to use the SLD style and elements to map the elevation value to the appropriate RGB value (with gradient interpolation between discrete values). For example:

          <ColorMap>
            <ColorMapEntry color="#000000" quantity="-10000" />
            <ColorMapEntry color="#FFFFFF" quantity="1667721.5" />
          </ColorMap>

It turns out that the above example does not span the full range of colors but rather creates gray values only. That is, it seems that it interpolate each color (red, green, blue) independent of each other.

Any idea how to make it interpolate values like that: #000000, #000001, #000002, ... , #0000FF, #000100, ..., #0001FF, ..., #FFFFFF?

Tx. [1]: https://docs.mapbox.com/data/tilesets/reference/mapbox-terrain-rgb-v1/

Answer 1

If you have the elevation data in the format you desire then that is the easiest option: it just works. However, if you want a more customized solution, here's what I've done for a project using the Mapbox Terrain-RGB format:

I have a scale of colors from dark blue to light blue to white. I want to be able to specify how many steps are used from light blue to white (default is 10). This code uses GDAL Python bindings. The following code snippet was used for testing. It just outputs the color mapping to a GeoTIFF file. To get values between 0 and 1, simply use value *= 1/num_steps. You can use that value in the lookup table to get an RGB value. If you're only interested in outputting the colors, you can ignore everything involving gdal_translate. The colors will automatically be stored in a single-band GeoTIFF. If you do want to re-use those colors, note that this version ignores alpha values (if present). You can use gdal_translate to add those. That code snippet is also available at my gist here.

import numpy as np
import gdal
from osgeo import gdal, osr

def get_color_map(num_steps):
    colors = np.zeros((num_steps, 3), dtype=np.uint8)
    colors[:, 0] = np.linspace(0, 255, num_steps, dtype=np.uint8)
    colors[:, 1] = colors[::-1, 0]
    return colors


ds = gdal.Open('/Users/myusername/Desktop/raster.tif')
band = ds.GetRasterBand(1) # Assuming single band raster
arr = band.ReadAsArray()
arr = arr.astype(np.float32)
arr *= 1/num_steps # Ensure values are between 0 and 1 (or use arr -= arr.min() / (arr.max() - arr.min()) to normalize to 0 to 1)


colors = get_color_map(num_steps) # Create color lookup table
colors[0] = [0, 0, 0] # Set black for no data so it doesn't show up as gray in the final product.


# Create new GeoTIFF with colors included (transparent alpha channel if possible). If you don't care about including the colors in the GeoTIFF, skip this.
cols = ds.RasterXSize
rows = ds.RasterYSize
out_ds = gdal.GetDriverByName('GTiff').Create('/Users/myusername/Desktop/raster_color.tif', cols, rows, 4)
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())
band = out_ds.GetRasterBand(1)
band.WriteArray(colors[arr]) # This can be removed if you don't care about including the colors in the GeoTIFF
band = out_ds.GetRasterBand(2)
band.WriteArray(colors[arr]) # This can be removed if you don't care about including the colors in the GeoTIFF
band = out_ds.GetRasterBand(3)
band.WriteArray(colors[arr]) # This can be removed if you don't care about including the colors in the GeoTIFF
band = out_ds.GetRasterBand(4)
alpha = np.zeros((rows, cols), dtype=np.uint8) # Create alpha channel to simulate transparency of no data pixels (assuming 0 is "no data" and non-zero is data). You can remove this if your elevation values are not 0.
alpha[arr == 0] = 255
band.WriteArray(alpha) # This can be removed if you don't care about including the colors in the GeoTIFF
out_ds.FlushCache()

This issue is also present in Rasterio when using a palette with multiple values. Here is an example.

However, if your raster has n-dimensions or is a masked array, the flip operation can be tricky. Here's a solution based on one of the answers in this stackoverflow question: How to vertically flip a 2D NumPy array?.

Answer 2

I'm trying to do the same with no luck, and i think it can't be done... Check this example. It's a "gradient" [-10000, -5000, -1000, -500 ... 100000000000000000, 5000000000000000000, 1000000000000000000] with the Mapbox color codification. The color progression/interpolation is anything but linear, so i think it can't be emulated in an SLD.