Graph illustrating the linear scaling and empirical quantile mapping (LSEQM) method for daily precipitation bias correction.

xkcd_lseqm.md

xkcd style for LSEQM illustration

This is example of Python script that could generate a xkcd graph-style to illustrate concept of Linear Scaling and Empirical Quantile Mapping (LSEQM) method for daily precipitation bias correction.

You can paste the below code into an online Python compiler like https://python-fiddle.com/ and grab the result instantly.

version a

version b

version c

flowchart

xkcd_lseqm_a.py

import matplotlib.pyplot as plt
import numpy as np

# Sample data (replace with your actual data)
obs_precip = [10, 5, 18, 12, 7, 15, 20]
sim_precip = [8, 3, 22, 15, 5, 12, 25]

# Function to calculate scaling factor
def scaling_factor(obs, sim):
    return sum(obs) / sum(sim)

# Linear scaling
scaling_factor_value = scaling_factor(obs_precip, sim_precip)
adjusted_precip = [p * scaling_factor_value for p in sim_precip]

# Combine data for CDF
data = sorted(obs_precip + adjusted_precip)

# Cumulative Distribution Function (CDF)
def cdf(data):
    cdf_values = []
    n = len(data)
    for i, d in enumerate(data):
        cdf_values.append((i + 1) / n)
    return cdf_values

obs_cdf = cdf(sorted(obs_precip))
adj_cdf = cdf(sorted(adjusted_precip))

# Bias correction using quantile mapping
def bias_correct(obs_val, obs_precip, adjusted_precip):
    obs_sorted = sorted(obs_precip)
    adj_sorted = sorted(adjusted_precip)
    obs_percentile = obs_sorted.index(obs_val) / len(obs_sorted)
    adj_index = int(obs_percentile * len(adj_sorted))
    return adj_sorted[adj_index]

# Sample day for correction (change index for different day)
day_to_correct = 5
obs_val = obs_precip[day_to_correct]
corr_precip = bias_correct(obs_val, obs_precip, adjusted_precip)

# Plotting in XKCD style
plt.xkcd()

plt.figure(figsize=(10, 6))

# Panel 1: Original Data
plt.subplot(221)
plt.bar(range(len(obs_precip)), obs_precip, color='blue', label='Observed')
plt.bar(range(len(sim_precip)), sim_precip, color='orange', label='Simulated')
plt.title('Original Precipitation')
plt.xlabel('Day')
plt.ylabel('Precipitation (mm)')
plt.legend()

# Panel 2: Linear Scaling
plt.subplot(222)
plt.bar(range(len(adjusted_precip)), adjusted_precip, color='green', label='Adjusted')
plt.title('Linear Scaling (a = {:.2f})'.format(scaling_factor_value))
plt.xlabel('Day')
plt.ylabel('Precipitation (mm)')
plt.legend()

# Panel 3: CDFs
plt.subplot(223)
plt.plot(sorted(obs_precip), np.linspace(0, 1, len(obs_precip), endpoint=False), marker='.', linestyle='none', label='Observed CDF', color='blue')
plt.plot(sorted(adjusted_precip), np.linspace(0, 1, len(adjusted_precip), endpoint=False), marker='.', linestyle='none', label='Adjusted CDF', color='red')
plt.title('Cumulative Distribution Functions (CDFs)')
plt.xlabel('Precipitation (mm)')
plt.ylabel('Probability')
plt.legend()

# Panel 4: Quantile Mapping (example)
plt.subplot(224)
plt.axvline(x=obs_precip[day_to_correct], color='red', linestyle='--', label='Observed ({})'.format(obs_precip[day_to_correct]))
plt.axvline(x=corr_precip, color='black', label='Bias-Corrected')
plt.axhline(y=obs_cdf[day_to_correct], color='red', linestyle='--')
plt.axhline(y=adj_cdf[int(obs_cdf[day_to_correct] * len(adj_cdf))], color='gray', linestyle='--')
plt.bar(day_to_correct, sim_precip[day_to_correct], color='orange', alpha=0.5, label='Original Psim')
plt.text(day_to_correct + 0.1, sim_precip[day_to_correct] + 0.5, '{:.1f}'.format(sim_precip[day_to_correct]), color='orange', ha='center', va='bottom')
plt.title('Quantile Mapping (Day {})'.format(day_to_correct))
plt.xlabel('Precipitation (mm)')
plt.ylabel('CDF Value')
plt.legend()

plt.tight_layout()
plt.show()

xkcd_lseqm_b.py

import matplotlib.pyplot as plt
import numpy as np

# Sample data (replace with your actual data)
obs_precip = [10, 5, 18, 12, 7, 15, 20]
sim_precip = [8, 3, 22, 15, 5, 12, 25]

# Function to calculate scaling factor
def scaling_factor(obs, sim):
    return sum(obs) / sum(sim)

# Linear scaling
scaling_factor_value = scaling_factor(obs_precip, sim_precip)
adjusted_precip = [p * scaling_factor_value for p in sim_precip]

# Combine data for CDF
data = sorted(obs_precip + adjusted_precip)

# Cumulative Distribution Function (CDF)
def cdf(data):
    cdf_values = []
    n = len(data)
    for i, d in enumerate(data):
        cdf_values.append((i + 1) / n)
    return cdf_values

obs_cdf = cdf(sorted(obs_precip))
adj_cdf = cdf(sorted(adjusted_precip))

# Bias correction using quantile mapping
def bias_correct(obs_val, obs_precip, adjusted_precip):
    obs_sorted = sorted(obs_precip)
    adj_sorted = sorted(adjusted_precip)
    obs_percentile = obs_sorted.index(obs_val) / len(obs_sorted)
    adj_index = int(obs_percentile * len(adj_sorted))
    return adj_sorted[adj_index]

# Correcting all days using quantile mapping
corr_precip_all = [bias_correct(val, obs_precip, adjusted_precip) for val in obs_precip]

# Sample day for correction (change index for different day)
day_to_correct = 5

# Plotting in XKCD style
plt.xkcd()

plt.figure(figsize=(15, 10))

# Panel 1: Original Data
plt.subplot(231)
width = 0.35
x = np.arange(len(obs_precip))
plt.bar(x - width/2, obs_precip, width, label='Observed', color='blue')
plt.bar(x + width/2, sim_precip, width, label='Simulated', color='orange')
plt.title('Original Precipitation Data')
plt.xlabel('Day')
plt.ylabel('Precipitation (mm)')
plt.legend()

# Panel 2: Mean Values
plt.subplot(232)
plt.bar(['Observed Mean', 'Simulated Mean'], [np.mean(obs_precip), np.mean(sim_precip)], color=['blue', 'orange'])
plt.title('Mean Values Before Scaling')
plt.ylabel('Precipitation (mm)')

# Panel 3: Linear Scaling
plt.subplot(233)
plt.bar(range(len(adjusted_precip)), adjusted_precip, color='green', label='Adjusted')
plt.title('Linear Scaling Applied')
plt.xlabel('Day')
plt.ylabel('Precipitation (mm)')
plt.legend()

# Panel 4: CDFs Comparison
plt.subplot(234)
plt.plot(sorted(obs_precip), np.linspace(0, 1, len(obs_precip), endpoint=False), marker='.', linestyle='none', label='Observed CDF', color='blue')
plt.plot(sorted(adjusted_precip), np.linspace(0, 1, len(adjusted_precip), endpoint=False), marker='.', linestyle='none', label='Adjusted CDF', color='red')
plt.title('Cumulative Distribution Functions (CDFs)')
plt.xlabel('Precipitation (mm)')
plt.ylabel('Probability')
plt.legend()

# Panel 5: Quantile Mapping (example)
plt.subplot(235)
plt.axvline(x=obs_precip[day_to_correct], color='red', linestyle='--', label='Observed ({})'.format(obs_precip[day_to_correct]))
plt.axvline(x=corr_precip_all[day_to_correct], color='black', label='Bias-Corrected')
plt.axhline(y=obs_cdf[day_to_correct], color='red', linestyle='--')
plt.axhline(y=adj_cdf[int(obs_cdf[day_to_correct] * len(adj_cdf))], color='gray', linestyle='--')
plt.bar(day_to_correct, sim_precip[day_to_correct], color='orange', alpha=0.5, label='Original Simulated')
plt.text(day_to_correct + 0.1, sim_precip[day_to_correct] + 0.5, '{:.1f}'.format(sim_precip[day_to_correct]), color='orange', ha='center', va='bottom')
plt.title('Quantile Mapping (Day {})'.format(day_to_correct))
plt.xlabel('Precipitation (mm)')
plt.ylabel('CDF Value')
plt.legend()

# Panel 6: Final Corrected Data
plt.subplot(236)
plt.bar(x - width/2, obs_precip, width, label='Observed', color='blue')
plt.bar(x + width/2, corr_precip_all, width, label='Bias-Corrected', color='red')
plt.title('Final Corrected Data')
plt.xlabel('Day')
plt.ylabel('Precipitation (mm)')
plt.legend()

plt.tight_layout()
plt.show()

xkcd_lseqm_flowchart.py

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# Function to add a box with text
def add_box(text, xy, ax, width=3, height=1, fontsize=18):
    box = mpatches.FancyBboxPatch(xy, width=width, height=height, boxstyle="round,pad=0.3", edgecolor='black', facecolor='white')
    ax.add_patch(box)
    ax.text(xy[0] + width / 2, xy[1] + height / 2, text, ha='center', va='center', fontsize=fontsize)

# Create the plot
plt.xkcd()
fig, ax = plt.subplots(figsize=(18, 8))

# Adding boxes for each step in the flowchart
steps = [
    "Set up\nPython env\n(xarray, numpy, etc)",
    "Load data as\nDataArray",
    "Linear Scaling",
    "Combine Data",
    "Calculate CDFs\nCumulative\nDistribution\nFunction",
    "Bias Correction\n(Quantile Mapping)",
    "Output\nDataArray",
    "Result"
]

positions = [
    (1, 6), (5, 6), (9, 6), (13, 6),
    (13, 3), (9, 3), (5, 3), (1, 3)
]

for step, pos in zip(steps, positions):
    add_box(step, pos, ax, fontsize=18)

# Adding arrows
arrows = [
    ((4, 6.5), (5, 6.5)),  # Step 1 to Step 2
    ((8, 6.5), (9, 6.5)),  # Step 2 to Step 3
    ((12, 6.5), (13, 6.5)),  # Step 3 to Step 4
    ((15, 6), (15, 4.2)),  # Step 4 to Step 5 (down)
    ((13, 3.5), (12, 3.5)),  # Step 5 to Step 6
    ((9, 3.5), (8, 3.5)),  # Step 6 to Step 7
    ((5, 3.5), (4, 3.5)),  # Step 7 to Step 8
]

for start, end in arrows:
    plt.arrow(start[0], start[1], end[0] - start[0], end[1] - start[1], head_width=0.2, head_length=0.2, fc='black', ec='black')

# Removing axes
ax.set_xlim(0, 18)
ax.set_ylim(2, 8)
ax.axis('off')

plt.tight_layout()
plt.show()

xkcd_lseqmdl.py

# Install required packages
!pip install -q requests matplotlib

import requests
import os
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
fm._load_fontmanager(try_read_cache=False)
from pathlib import Path
import warnings
import numpy as np

# Set the figure DPI to 300
plt.rcParams['figure.dpi'] = 300

# Suppress font-related warnings
warnings.filterwarnings("ignore", 
                       message="findfont: Font family .* not found", 
                       category=UserWarning, 
                       module="matplotlib.font_manager")

# Download the xkcd font file
font_url = "https://github.com/ipython/xkcd-font/raw/master/xkcd-script/font/xkcd-script.ttf"
font_path = "/tmp/xkcd-script.ttf"

# Download the font
response = requests.get(font_url)
with open(font_path, "wb") as f:
    f.write(response.content)

# Update matplotlib's font cache
fm.fontManager.addfont(font_path)

# Use XKCD style for a fun illustration
plt.xkcd()

# Increase the global font size
plt.rcParams.update({
    'font.size': 16,              # Increase base font size
    'axes.titlesize': 18,         # Title font size
    'axes.labelsize': 14,         # Axis label font size
    'xtick.labelsize': 14,        # X-tick label size
    'ytick.labelsize': 14,        # Y-tick label size
    'legend.fontsize': 16,        # Legend font size
    'figure.titlesize': 24        # Figure title size
})

# Sample data (replace with your actual data)
obs_precip = [10, 5, 18, 12, 7, 15, 20]
sim_precip = [8, 3, 22, 15, 5, 12, 25]

# ----------------------------
# Step 1: Linear Scaling (LS)
# ----------------------------
def scaling_factor(obs, sim):
    return sum(obs) / sum(sim)

ls_factor = scaling_factor(obs_precip, sim_precip)
# LS adjusted precipitation is the simulated data scaled by LS factor
adjusted_precip = [p * ls_factor for p in sim_precip]

# ----------------------------
# Step 2: Quantile Mapping (EQM)
# ----------------------------
# Compute cumulative distribution function (CDF)
def cdf(data):
    n = len(data)
    return [(i + 1) / n for i in range(n)]

# Compute CDFs for observed and LS-adjusted data
obs_sorted = sorted(obs_precip)
adjusted_sorted = sorted(adjusted_precip)
obs_cdf = cdf(obs_sorted)
adj_cdf = cdf(adjusted_sorted)

# EQM: bias correction via quantile mapping
def bias_correct(obs_val, obs_precip, adjusted_precip):
    obs_sorted = sorted(obs_precip)
    adj_sorted = sorted(adjusted_precip)
    obs_percentile = obs_sorted.index(obs_val) / len(obs_sorted)
    # Ensure index is within bounds
    adj_index = min(int(obs_percentile * len(adj_sorted)), len(adj_sorted) - 1)
    return adj_sorted[adj_index]

# Compute bias-corrected precipitation using EQM
corr_precip_all = [bias_correct(val, obs_precip, adjusted_precip) for val in obs_precip]

# ----------------------------
# Step 3: Deep Learning (DL) Correction
# ----------------------------
# Simulated DL adjustment: for extreme events (e.g., precipitation > 15), apply an additional factor.
def apply_deep_learning_adjustment(precip_data, threshold=15, adjustment_factor=1.1):
    corrected_data = []
    for val in precip_data:
        if val > threshold:
            corrected_data.append(val * adjustment_factor)
        else:
            corrected_data.append(val)
    return corrected_data

# Apply DL correction on the EQM corrected values
dl_corrected_precip_all = apply_deep_learning_adjustment(corr_precip_all, threshold=15, adjustment_factor=1.1)

# Compute DL adjustment differences: DL correction minus EQM correction
dl_diff = [dl - eqm for dl, eqm in zip(dl_corrected_precip_all, corr_precip_all)]

# ----------------------------
# Step 4: Define a Sample Day for Detailed Comparison
# ----------------------------
day_to_correct = 5  # change index if needed

# ----------------------------
# Step 5: Composite Plot: 2 Rows × 4 Columns (Figures 1 to 8)
# ----------------------------

# Increase figure size to give more room for text
fig, axes = plt.subplots(2, 4, figsize=(24, 12))

# --- Figure 1: Original Precipitation Data ---
ax = axes[0, 0]
width = 0.35
x = np.arange(len(obs_precip))
ax.bar(x - width/2, obs_precip, width, label='Observed', color='blue')
ax.bar(x + width/2, sim_precip, width, label='Satellite', color='orange')
ax.set_title('1: Precipitation Data')
ax.set_xlabel('Day')
ax.set_ylabel('Precip (mm)')
ax.legend()

# --- Figure 2: LS – Mean Values Before Scaling ---
ax = axes[0, 1]
ax.bar(['Observed Mean', 'Satellite Mean'], 
       [np.mean(obs_precip), np.mean(sim_precip)], 
       color=['blue', 'orange'])
ax.set_title('2: Mean Values Before Scaling (LS)')
ax.set_ylabel('Precip (mm)')

# --- Figure 3: LS – Linear Scaling Applied ---
ax = axes[0, 2]
ax.bar(range(len(adjusted_precip)), adjusted_precip, color='green', label='LS Adjusted')
ax.set_title('3: Linear Scaling Applied')
ax.set_xlabel('Day')
ax.set_ylabel('Precip (mm)')
ax.legend()

# --- Figure 4: LS/EQM – CDFs Comparison ---
ax = axes[0, 3]
# Plot observed and LS-adjusted sorted values vs. their CDF
ax.plot(obs_sorted, np.linspace(0, 1, len(obs_sorted), endpoint=False), 
        marker='.', markersize=12, linestyle='none', label='Observed CDF', color='blue')
ax.plot(adjusted_sorted, np.linspace(0, 1, len(adjusted_sorted), endpoint=False),
        marker='.', markersize=12, linestyle='none', label='LS Adjusted CDF', color='red')
ax.set_title('4: CDFs Comparison (LS/EQM)')
ax.set_xlabel('Precip (mm)')
ax.set_ylabel('Probability')
ax.legend()

# --- Figure 5: EQM – Quantile Mapping Example ---
ax = axes[1, 0]
# Mark the observed value and the EQM-corrected value for the sample day
ax.axvline(x=obs_precip[day_to_correct], color='red', linestyle='--', linewidth=2,
           label=f'Observed ({obs_precip[day_to_correct]})')
ax.axvline(x=corr_precip_all[day_to_correct], color='black', linewidth=2,
           label='Bias-Corrected (EQM)')
ax.axhline(y=obs_cdf[day_to_correct], color='red', linestyle='--', linewidth=2)
# Plot original simulated value for reference
ax.bar(day_to_correct, sim_precip[day_to_correct], color='orange', alpha=0.5, 
       label='Original Satellite')
ax.text(day_to_correct + 0.1, sim_precip[day_to_correct] + 0.005, 
        f'{sim_precip[day_to_correct]:.1f}', color='orange', fontsize=14, ha='center', va='bottom')
ax.set_title(f'5: Quantile Mapping Example (Day {day_to_correct})')
ax.set_xlabel('Precip (mm)')
ax.set_ylabel('CDF Value')
ax.legend()

# --- Figure 6: DL – DL Adjustment Differences ---
ax = axes[1, 1]
ax.bar(range(len(dl_diff)), dl_diff, color='purple')
ax.set_title('6: DL Adjustment Differences\n(DL Correction minus EQM)')
ax.set_xlabel('Day')
ax.set_ylabel('Difference (mm)')

# --- Figure 7: DL – Scatter Plot of EQM vs. DL Corrected Values ---
ax = axes[1, 2]
ax.scatter(corr_precip_all, dl_corrected_precip_all, color='magenta', s=100)  # Increased marker size
ax.plot([min(corr_precip_all), max(corr_precip_all)], 
        [min(corr_precip_all), max(corr_precip_all)], 
        color='gray', linestyle='--', linewidth=2)  # thicker diagonal reference line
ax.set_title('7: EQM vs. DL Corrected Values')
ax.set_xlabel('Bias-Corrected (EQM)')
ax.set_ylabel('DL Corrected')

# --- Figure 8: Final Corrected Data with DL ---
ax = axes[1, 3]
ax.bar(x - width/2, obs_precip, width, label='Observed', color='blue')
ax.bar(x + width/2, dl_corrected_precip_all, width, label='Bias-Corrected + DL', color='red')
ax.set_title('8: Final Corrected Data')
ax.set_xlabel('Day')
ax.set_ylabel('Precip (mm)')
ax.legend()

# Add more space between subplots to prevent text overlap
plt.subplots_adjust(wspace=0.3, hspace=0.4)

plt.tight_layout()
plt.show()