I’ve written some code to create subplots with maps based on the JRC Global Surface Water dataset. My code creates the maps just fine, but when I want to add a legend to explain the colors, they look off. I’ve written extra code to test the colorbars and the legends, and for some reason, the colors never fill the bar fully.

I assume this is due to normalization, but when I don’t normalize the categorical subplots, they don’t plot at all.

import pystac_client
import rasterio
import pystac
import planetary_computer
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize, ListedColormap
import logging
import dask.array as da
from dask.distributed import Client, LocalCluster
import rioxarray
import xarray as xr
import rasterio
from rasterio.enums import Resampling

# Set up logging
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

city_name="Manaus"
# Set up Dask client with memory limit
cluster = LocalCluster(memory_limit='5GB')
client = Client(cluster)
# Set up the STAC client
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

bbox_of_interest = (-62.9825041, -6.1316333, -56.9825041, -0.13163329999999984)
area_of_interest = {
    "type": "Polygon",
    "coordinates": [[
        [-56.9825041, -0.13163329999999984],
        [-56.9825041, -6.1316333],
        [-62.9825041, -6.1316333],
        [-62.9825041, -0.13163329999999984],
        [-56.9825041, -0.13163329999999984]
    ]]
}

# Search for items
search = catalog.search(collections=["jrc-gsw"], intersects=area_of_interest)
items = list(search.get_items())  # Retrieve all items as a list

print(f"Returned {len(items)} items")
item = items[0]

print(f"Processing Item ID: {item.id}")
cog_assets = [
    asset_key for asset_key, asset in item.assets.items()
    if asset.media_type == pystac.MediaType.COG
]

ASSET_RANGES = {
    'change': (0, 200),
    'extent': (0, 1),
    'occurrence': (0, 100),
    'recurrence': (0, 100),
    'seasonality': (1, 12),
    'transitions': (0, 10)
}

# Create colormaps for each asset
cmaps = {}
for asset_key in cog_assets:
    with rasterio.open(item.assets[asset_key].href) as src:
        colormap_def = src.colormap(1)
        colormap = [np.array(colormap_def[i]) / 256 for i in range(256)]
cmaps[asset_key] = ListedColormap(colormap)

def analyze_asset(item, asset_key):
    try:
        with rasterio.open(item.assets[asset_key].href) as src:
            data = src.read(1)
            colormap_def = src.colormap(1)
            colormap = [np.array(colormap_def[i]) / 255 for i in range(256)]
            
            # Use the predefined range for this asset type
            min_value, max_value = ASSET_RANGES.get(asset_key, (data.min(), data.max()))
            
            return {
                'colormap': colormap,
                'min_value': min_value,
                'max_value': max_value
            }
    except Exception as e:
        logger.error(f"Error analyzing asset {asset_key}: {str(e)}")
        return None

# Analyze each asset
results = {}
for asset_key in cog_assets:
    result = analyze_asset(items[0], asset_key)
    if result:
        results[asset_key] = result
        logger.info(f"nAsset: {asset_key}")

fig, axs = plt.subplots(len(results), 1, figsize=(10, 2*len(results)))
for i, (asset_key, result) in enumerate(results.items()):
    cmap = ListedColormap(result['colormap'])
    min_val, max_val = result['min_value'], result['max_value']
    
    # Create an array of values from min to max
    data_range = np.linspace(min_val, max_val, 256)
    
    # Plot the range of values for this asset type
    im = axs[i].imshow(data_range.reshape(1, -1), cmap=cmap, aspect='auto', extent=[min_val, max_val, 0, 1])
    axs[i].set_title(f"{asset_key} colormap")
    axs[i].set_yticks([])
    
    # Add colorbar
    cbar = plt.colorbar(im, ax=axs[i], orientation='horizontal', aspect=25)
    
    # Set appropriate tick locations
    if asset_key == 'extent':
        cbar.set_ticks([0, 1])
    elif asset_key == 'seasonality':
        cbar.set_ticks(range(1, 13))
    elif asset_key == 'transitions':
        cbar.set_ticks(range(0, 11))
    else:
        cbar.set_ticks(np.linspace(min_val, max_val, 5))

plt.tight_layout()
plt.show()

I’ve checked whether it’s due to normalizing, tried mapping the values, tried setting ranges.