This is the code to make the figure tweeted in https://twitter.com/rrealrangel/status/1336797063724617728. It is product of personal work and originally was not intended to be published, so it is not as clear as it could be with a little more time available. Thanks for your interest. ;)

cmap_figure

# -*- coding: utf-8 -*-
"""
Created on Thu Dec 10 16:49:05 2020
@author: rreal
"""
import matplotlib as mpl
import numpy as np
import pandas as pd

BDCN_TS = "C://..."   # Full path of observed data CSV input files.
CHIRPS_TS_RAW = "C://..."   # Full path of simulated data CSV input files.
CHIRPS_TS_COR = "C://..."   # Full path of corrected data CSV input files.

def compare_values(output_file):
    def mae(mod, obs):
        """Compute the mean absolute error (MAE)."""
        return (mod - obs).abs().mean()

    def mbe(mod, obs):
        """Compute the mean bias error (MBE)."""
        return (mod - obs).mean()
    cor = pd.read_csv(
        filepath_or_buffer=CHIRPS_TS_COR, index_col=0, parse_dates=True
        ).stack()
    obs = pd.read_csv(
        filepath_or_buffer=BDCN_TS, index_col=0, parse_dates=True
        ).stack().loc[cor.index]
    mod = pd.read_csv(
        filepath_or_buffer=CHIRPS_TS_RAW, index_col=0, parse_dates=True
        ).stack().loc[cor.index]

    # 2D histograms: configuration.
    xmax = ymax = 830
    bin_size = 15

    # 2D histograms: observed data versus modeled data.
    obs_mod_hist2d, xedges, yedges = np.histogram2d(
        x=obs,
        y=mod,
        bins=np.arange(0, xmax, bin_size)
        )
    obs_mod_ltrend_f = np.poly1d(np.polyfit(obs, mod, 1))
    obs_mod_ltrend = obs_mod_ltrend_f(obs)

    # 2D histograms: observed data versus corrected modeled data.
    obs_cor_hist2d, xedges, yedges = np.histogram2d(
        x=obs,
        y=cor,
        bins=np.arange(0, xmax, bin_size)
        )
    obs_cor_ltrend_f = np.poly1d(np.polyfit(obs, cor, 1))
    obs_cor_ltrend = obs_cor_ltrend_f(obs)

    extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

    # Define a new color map.
    cmap = mpl.pyplot.cm.inferno_r
    cmaplist = [cmap(i) for i in range(cmap.N)]
    cmaplist[0] = (1, 1, 1, 1)  # force 1st color entry to be white.

    for i in range(1, 50):  # force a few 2nd color entries to be light grey.
        cmaplist[i] = (0.85, 0.85, 0.85, 1)

    cmap = mpl.colors.LinearSegmentedColormap.from_list(
        'Custom cmap', cmaplist, cmap.N
        )
    bounds = np.sort(np.concatenate((
        np.array([0.99999999999, 5]), np.linspace(0, 100, 6)
        )))
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

    # Plot the data.
    mpl.pyplot.rcParams["figure.figsize"] = (3.937, 5.906)
    fig, ax = mpl.pyplot.subplots(nrows=1, ncols=2, sharey=True)

    # Left panel.
    ax[0].plot([0, xmax],   [0, ymax],      c='k',  lw=0.5)
    ax[0].plot(obs, obs_mod_ltrend, c='r',  lw=0.5)
    hm = ax[0].imshow(
        X=obs_mod_hist2d.T,
        extent=extent,
        cmap=cmap,
        norm=norm,
        origin='lower',
        interpolation="nearest"
        )

    # Right panel.
    ax[1].plot([0, xmax],   [0, ymax],          c='k',  lw=0.5)
    ax[1].plot(obs, obs_cor_ltrend, c='r',  lw=0.5)
    ax[1].imshow(
        X=obs_cor_hist2d.T,
        extent=extent,
        cmap=cmap,
        norm=norm,
        origin='lower',
        interpolation="nearest"
        )

    # Axes labels.
    font = {'family': 'Source Sans Pro', 'weight': 'normal', 'size': 9}
    ax[0].set_xlim(right=xmax)
    ax[0].set_ylim(top=ymax)
    ax[0].set_xlabel(xlabel='Precip. observada (mm)', fontdict=font)
    ax[1].set_xlabel(xlabel='Precip. observada (mm)', fontdict=font)
    ax[0].set_ylabel(ylabel='Precip. CHIRPS (mm)', fontdict=font)

    for axis in ax:
        for tick in axis.get_xticklabels():
            tick.set_fontname('Source Sans Pro')
            tick.set_fontsize(9)

        for tick in axis.get_yticklabels():
            tick.set_fontname('Source Sans Pro')
            tick.set_fontsize(9)

    # Annotations.
    mae_1 = mae(mod=mod, obs=obs)
    mbe_1 = mbe(mod=mod, obs=obs)
    text = "MAE = {:.3f}\nMBE = {:.3f}".format(mae_1, mbe_1)
    ax[0].text(
        x=0.05, y=0.80, s=text, fontdict=font, transform=ax[0].transAxes
        )
    mae_2 = mae(mod=cor, obs=obs)
    mbe_2 = mbe(mod=cor, obs=obs)
    text = "MAE = {:.3f}\nMBE = {:.3f}".format(mae_2, mbe_2)
    ax[1].text(
        x=0.05, y=0.80, s=text, fontdict=font, transform=ax[1].transAxes
        )

    # Colorbar.
    cbar = fig.colorbar(hm, ax=ax, location='bottom', spacing='proportional')
    cbar_ticks = np.sort(
        np.concatenate([[1, 5], np.linspace(0, 100, 6)[1:-1]])
        ).astype(int)
    cbar.set_ticks(cbar_ticks)
    cbar.set_ticklabels(cbar_ticks)
    cbar.set_label('Frecuencia', size=9, family='Source Sans Pro')

    for tick in cbar.ax.get_xticklabels():
        tick.set_fontname('Source Sans Pro')
        tick.set_fontsize(9)

    # Export the figure.
    mpl.pyplot.tight_layout()
    fig.savefig(fname=output_file, dpi=400, bbox_inches='tight')