Can add_raster work with a matplotlib.colors.LinearSegmentedColormap?

[nve-sak]

Has anyone succeeded in feeding a matplotlib.colors.LinearSegmentedColormap to add_raster(), or know that it's not supported?

I need to display a raster image of continuous data (in a numpy array) as reclassified into danger levels green-yellow-orange-red.

The colours vary largely in value interval size, so an even-spaced matplotlib.colors.ListedColormap won't do. Unfortunately the two types have rather different formulation, and add_raster fails to produce the colours I expect. There is no crash or error.

Of course I could resample the numpy array, but I am letting the user browse through a long time series of maps, so would prefer not.

[giswqs]

Try updating Localtileserver to the latest version. It should work. See banesullivan/localtileserver#223 (comment)

[nve-sak]

Thanks for the quick answer :-) I had localtileserver 0.10.4 dated October 4, upgraded to 0.10.5 dated October 27, to no effect. But now I know it's intended to work, I'll cut a bare-bone version to test :-)

[nve-sak]

Code to reconstruct the unexpected behaviour without relying on any files:
`
import os, numpy as np
from datetime import datetime
from leafmap import Map, array_to_memory_file
from pprint import pprint
from matplotlib.colors import LinearSegmentedColormap
del os.environ['PROJ_LIB'] # Too many versions of proj hanging around...

Create the raster to plot:

l = list(range(0,80,1)) # Values to be classified into green, yellow, orange, red danger level
l.extend([65535]*10) # 65535 means missing data
l.extend([65533]*10) # 65533 is another special value
nparr = np.array(l, dtype=np.float32)
nparr[nparr<65000] /= 10
nparr = nparr.reshape([10,10])
GT33 = (1000.0,0.0,300000.0,0.0,-1000.0,7000000.0)
rstimg = array_to_memory_file(array=nparr, cellsize=1000, crs="epsg:32633", transform=GT33) # UTM zone 33N

Create a colormap according to https://matplotlib.org/stable/api/_as_gen/matplotlib.colors.LinearSegmentedColormap.html

valscale = 65535.0 # Values are on a scale from 0 to 65535, we need to normalise.
colscale = 255.0 # This dict is hard-coded here, in real life it is imported from a file I cannot alter.
cmapdict = {
'blue': [[0.0, 96/colscale, 96/colscale], # Plot as green
[2.0/valscale, 96/colscale, 0.0 ], # Plot as yellow
[4.0/valscale, 0.0, 0.0 ], # Plot as orange
[6.0/valscale, 0.0, 0.0 ], # Plot as red
[65533/valscale, 0.0, 229/colscale], # Plot as grey (special value)
[1.0, 229/colscale, 229/colscale] ] # Plot as transparent (missing data)
'green': [[0.0, 212/colscale, 212/colscale],
[2.0/valscale, 212/colscale, 252/colscale],
[4.0/valscale, 252/colscale, 196/colscale],
[6.0/valscale, 196/colscale, 0.0 ],
[65533/valscale, 0.0, 229/colscale],
[1.0, 229/colscale, 229/colscale] ],
'red': [ [0.0, 40/colscale, 40/colscale],
[2.0/valscale, 40/colscale, 248/colscale],
[4.0/valscale, 248/colscale, 248/colscale],
[6.0/valscale, 248/colscale, 248/colscale],
[65533/valscale, 248/colscale, 229/colscale],
[1.0, 229/colscale, 229/colscale] ] }
cmap = LinearSegmentedColormap("mycolordict", cmapdict)
pprint(cmap)
pprint(cmap.dict)

M = Map(basemap='OpenStreetMap', center=[65, 18], zoom=5, scroll_wheel_zoom=True)
M.add_raster(rstimg, colormap=cmap, layer_name = "Test", nodata=65535, opacity=0.9)
M
`

This result in the following image:

The main part of the shown rectangle is 8 green rows, should have been 2 green, 2 yelllow, 2 orange, 2 red.
Transparent row 9 is correct, these are missing values. Red row 10 should have been grey.

I am trying to mimic something like this (from www.xgeo.no):
https://www.xgeo.no/?p=flomogjordskred&m=bmGeodataGraatone%3BMapLayer_gwb_qttprqttxyrx30yr&l=en&d=1691438400000&e=-1226514%7C6369094%7C2230602%7C8002826

[giswqs]

Or try add_cog_layer if the raster is available through http

https://geog-312.gishub.org/book/geospatial/leafmap.html#using-a-custom-colormap

[nve-sak]

Thanks again. I extract the numpy arrays from local files which I need to manipulate a bit, so I hope I can stick to numpy arrays :-)

[banesullivan]

You're close! What I recommend is starting simple with a plot in something like matplotlib. This will help you properly create the colormap without navigating other bugs or oddities of the more complicated geospatial stack like localtilserver.

This is what I have:

import os, numpy as np
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt

# Create the raster to plot:
l = list(range(0,80,1)) # Values to be classified into green, yellow, orange, red danger level
l.extend([65535]*10) # 65535 means missing data
l.extend([65533]*10) # 65533 is another special value
nparr = np.array(l, dtype=np.float32)
nparr[nparr<65000] /= 10
nparr = nparr.reshape([10,10])
GT33 = (1000.0,0.0,300000.0,0.0,-1000.0,7000000.0)
nparr

# So what we have here are valid clim ~0-10, 65535 is NaN, and otherwise outside of bounds = grey

colors = ['green', 'yellow', 'orange', 'red']
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', colors)
cmap.set_under('grey')
cmap.set_over('grey')
cmap.set_bad((0, 0, 0, 0))  # set no data as transparent
im = nparr.copy()
im[np.where(im == 65535)] = np.nan  # this is what `nodata` does in localtileserver
plt.imshow(im, cmap=cmap, clim=(0,10))
plt.colorbar()

Now we have a colormap built properly. We should be able to pass this directly to localtileserver for geospatial visualization too. However, this colormap isn't serializing nicely for localtilerver and its unable to display. I'm going to open an issue in localtileserver for this and see what I can do to resolve it.

[banesullivan]

See banesullivan/localtileserver#231

[nve-sak]

Thanks a lot for diving into this. In pure Matplotlib, I have also managed to get what I want in two ways:

1. Using matplotlib.colors.from_levels_and_colors,
2. Using matplotlib.colors.LinearSegmentedColormap and adding your set_bad and 65535_to_nan lines.

Note that I am after a classifiied color map, not a gradual green-to-red transition.

But as you say, it seems that LeafMap does not read the color maps well, my suspicion is that the normalisation info is not picked up. Try removing the 'N=65535" from method 1 and the norm=norm3 from method 2.

As for the question, I consider it answered, thanks a lot :-)

Here is what i have after your suggestions:

valscale = 65535.0         # Values are on a scale from 0 to 65535, we need to normalise.
colscale = 255.0           # This dict is hard-coded here, in real life it is imported from one of
cmapdict = {               # hundreds of similar files for differernt maps, which I cannot modify.
 'red':   [[0.0,             40/colscale,  40/colscale],    # [threshold, left RGB component, right RGB component]
           [2.0/valscale,    40/colscale, 248/colscale],    # Left in this row matches right in the row above (lower category)
           [4.0/valscale,   248/colscale, 248/colscale],    # to have a constant color spanning an interval. 
           [6.0/valscale,   248/colscale, 248/colscale],    # Left != Right means that this threshold breaks color.
           [65533/valscale, 248/colscale, 229/colscale],
           [1.0,            229/colscale, 1.0] ], 
 'green': [[0.0,            212/colscale, 212/colscale],
           [2.0/valscale,   212/colscale, 252/colscale],
           [4.0/valscale,   252/colscale, 196/colscale],
           [6.0/valscale,   196/colscale, 0.0         ],
           [65533/valscale, 0.0,          229/colscale],
           [1.0,            229/colscale, 1.0] ],
 'blue':  [[0.0,             96/colscale,  96/colscale],    # green  [0.0 to 2.0>     (RGB = 40 212 96)
           [2.0/valscale,    96/colscale, 0.0         ],    # yellow [2.0 to 4.0>     (RGB = 248 252 0)
           [4.0/valscale,   0.0,          0.0         ],    # orange [4.0 to 6.0>     (RGB = 248 196 0)
           [6.0/valscale,   0.0,          0.0         ],    # red    [6.0 to 65533>   (RGB = 248 0 0)
           [65533/valscale, 0.0,          229/colscale],    # grey   [65533 to 65535> (RGB = 229 229 229)
           [1.0,            229/colscale, 1.0] ] } #  missing [65535>   (transparent)
cmap = mcolors.LinearSegmentedColormap("mycolordict", segmentdata=cmapdict, N=65535)
cmap.set_bad((0, 0, 0, 0))  # set no data as transparent
im = nparr.copy()
im[np.where(im == 65535)] = np.nan  # this is what `nodata` does in localtileserver
print(type(cmap))
plt.imshow(im, cmap=cmap)  # This gives the desired result

cs = 255.0
levels = [0, 2, 4, 6, 65533, 65535]
colors = [(40./cs,212./cs,96./cs), (248./cs,252./cs,0./cs), (248./cs,196./cs,0./cs), 
          (248./cs, 0./cs,0./cs), (229./cs,229./cs,229./cs)]
cmap3, norm3 = mcolors.from_levels_and_colors(levels, colors)
print(type(cmap3))
print(type(norm3))
plt.imshow(nparr, cmap=cmap3, norm=norm3)  # This gives the desired result

from leafmap import Map, array_to_image
if 'PROJ_LIB' in os.environ:
    del os.environ['PROJ_LIB']      # Too many versions of proj hanging around...
GT33 = (1000.0,0.0,300000.0,0.0,-1000.0,7000000.0)  # 1*1 km resolution, somewhere southeast of Trondheim, Norway
rstimg = array_to_image(array=nparr, cellsize=1000, crs="epsg:32633", transform=GT33) # UTM zone 33N
cur_center = [63, 11]
cur_zoom = 8
lMap = Map(basemap='OpenStreetMap', center=cur_center, zoom=cur_zoom, scroll_wheel_zoom=True)
lMap.add_raster(rstimg, colormap=cmap, layer_name = "Test", nodata=1., opacity=0.9, vmin=0.0, vmax=0.999999)#, norm=norm4)
lMap     # The image is turned upside down, that's fixable. The colorscale is the problem.

lMap2 = Map(basemap='OpenStreetMap', center=cur_center, zoom=cur_zoom, scroll_wheel_zoom=True)
lMap2.add_raster(rstimg, colormap=cmap3, layer_name = "Test", nodata=1., opacity=0.9, vmin=0.0, vmax=0.999999)#, norm=norm4)
lMap2   # Again the upside down image is fixable. The colorscale is different, but still not correct.