Change Colorbar Scaling in Matplotlib

Question

Using NASA's SRTM data, I've generated a global elevation heatmap.

The problem is, however, the continents tend to blend in with the ocean because of the range of elevation values. Is it possible to change the colorbar's scale so that the edges of the continents are more distinct from the ocean? I've tried different cmaps, but they all seem to suffer from the problem.

Here is my code. I'm initializing a giant array (with 0s) to hold global elevation data, and then populating it file by file from the SRTM dataset. Each file is 1 degree latitude by 1 degree longitude.

Another question I had was regarding the map itself. For some reason, the Appalachian Mountains seem to have disappeared entirely.

import os
import numpy as np
from .srtm_map import MapGenerator
from ..utils.hgt_parser import HGTParser
from tqdm import tqdm
import cv2
import matplotlib.pyplot as plt
import richdem as rd

class GlobalMapGenerator():
    def __init__(self):
        self.gen = MapGenerator()
        self.base_dir = "data/elevation/"
        self.hgt_files = os.listdir(self.base_dir)
        self.global_elevation_data = None
    
    def shrink(data, rows, cols):
        return data.reshape(rows, data.shape[0]/rows, cols, data.shape[1]/cols).sum(axis=1).sum(axis=2)

    def GenerateGlobalElevationMap(self, stride):
        res = 1201//stride
        max_N = 59
        max_W = 180
        max_S = 56
        max_E = 179
        
        # N59 --> N00
        # S01 --> S56
        # E000 --> E179
        # W180 --> W001
        
        # Initialize array global elevation
        self.global_elevation_data = np.zeros(( res*(max_S+max_N+1), res*(max_E+max_W+1) ))

        print("Output Image Shape:", self.global_elevation_data.shape)

        for hgt_file in tqdm(self.hgt_files):
            lat_letter = hgt_file[0]
            lon_letter = hgt_file[3]
            lat = int(hgt_file[1:3])
            lon = int(hgt_file[4:7])

            if lat_letter == "S":
                # Shift south down by max_N, but south starts at S01 so we translate up by 1 too
                lat_trans = max_N + lat - 1
            else:
                # Bigger N lat means further up. E.g. N59 is at index 0 and is higher than N00
                lat_trans = max_N - lat
            
            if lon_letter == "E":
                # Shift east right by max_W
                lon_trans = max_W + lon
            else:
                # Bigger W lon means further left. E.g. W180 is at index 0 and is more left than W001
                lon_trans = max_W - lon

            # load in data from file as resized
            data = cv2.resize(HGTParser(os.path.join(self.base_dir, hgt_file)), (res, res))
            
            # generate bounds (x/y --> lon.lat for data from this file for the giant array)
            lat_bounds = [res*lat_trans, res*(lat_trans+1)]
            lon_bounds = [res*lon_trans, res*(lon_trans+1)]
            
            try:
                self.global_elevation_data[ lat_bounds[0]:lat_bounds[1],  lon_bounds[0]:lon_bounds[1] ] = data
            except:
                print("REFERENCE ERROR: " + hgt_file)
                print("lat: ", lat_bounds)
                print("lon: ", lon_bounds)

        # generate figure
        plt.figure(figsize=(20,20))
        plt.imshow(self.global_elevation_data, cmap="rainbow")
        plt.title("Global Elevation Heatmap")
        plt.colorbar()
        plt.show()
        np.save("figures/GlobalElevationMap.npy", self.global_elevation_data)
        plt.savefig("figures/GlobalElevationMap.png")
    
    def GenerateGlobalSlopeMap(self, stride):
        pass

Answer 1

Use a TwoSlopeNorm (docs) for your norm, like the example here.

From the example:

Sometimes we want to have a different colormap on either side of a conceptual center point, and we want those two colormaps to have different linear scales. An example is a topographic map where the land and ocean have a center at zero, but land typically has a greater elevation range than the water has depth range, and they are often represented by a different colormap.

If you set the midpoint at sea level (0), then you can have two very different scalings based on ocean elevation vs land elevation.

Example code (taken from the example linked above):

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cbook as cbook
from matplotlib import cm

dem = cbook.get_sample_data('topobathy.npz', np_load=True)
topo = dem['topo']
longitude = dem['longitude']
latitude = dem['latitude']

fig, ax = plt.subplots()
# make a colormap that has land and ocean clearly delineated and of the
# same length (256 + 256)
colors_undersea = plt.cm.terrain(np.linspace(0, 0.17, 256))
colors_land = plt.cm.terrain(np.linspace(0.25, 1, 256))
all_colors = np.vstack((colors_undersea, colors_land))
terrain_map = colors.LinearSegmentedColormap.from_list(
    'terrain_map', all_colors)

# make the norm:  Note the center is offset so that the land has more
# dynamic range:
divnorm = colors.TwoSlopeNorm(vmin=-500., vcenter=0, vmax=4000)

pcm = ax.pcolormesh(longitude, latitude, topo, rasterized=True, norm=divnorm,
                    cmap=terrain_map, shading='auto')
# Simple geographic plot, set aspect ratio beecause distance between lines of
# longitude depends on latitude.
ax.set_aspect(1 / np.cos(np.deg2rad(49)))
ax.set_title('TwoSlopeNorm(x)')
cb = fig.colorbar(pcm, shrink=0.6)
cb.set_ticks([-500, 0, 1000, 2000, 3000, 4000])
plt.show()

See how it scales numbers with this simple usage (from docs):

>>> import matplotlib. Colors as mcolors
>>> offset = mcolors.TwoSlopeNorm(vmin=-4000., vcenter=0., vmax=10000)
>>> data = [-4000., -2000., 0., 2500., 5000., 7500., 10000.]
>>> offset(data)
array([0., 0.25, 0.5, 0.625, 0.75, 0.875, 1.0])