Joyplot of dem - Have fun

joy.py

import numpy as np
import matplotlib.pyplot as plt
import rasterio
from scipy.ndimage import gaussian_filter


input_dem="input.tif"
skip = 250
line_spacing = 900
vertical_exaggeration = 3.5
pad_width = 500
baseline = 10
edge_jitter = 25
body_jitter = 0.3
fade_len = 250 


def taper_from_edge_value(y, pad_width, baseline=0, edge_jitter=1.0, body_jitter=0.3, fade_len=50):
    def random_walk_noise(length, jitter_strength=1.0):
        walk = np.cumsum(np.random.normal(scale=jitter_strength, size=length))
        walk -= np.mean(walk)
        return walk

    left_val = y[0]
    right_val = y[-1]

    ramp = 0.5 * (1 - np.cos(np.linspace(0, np.pi, pad_width)))

    
    left = baseline + (left_val - baseline) * ramp
    right = baseline + (right_val - baseline) * ramp[::-1]

    
    left += random_walk_noise(pad_width, jitter_strength=edge_jitter)
    right += random_walk_noise(pad_width, jitter_strength=edge_jitter)

    
    y_center = y + np.random.normal(0, body_jitter, size=y.shape)

    
    weights = np.linspace(0, 1, fade_len)

    y_center[:fade_len] = (
        weights * y_center[:fade_len] +
        (1 - weights) * left[-fade_len:]
    )
    y_center[-fade_len:] = (
        weights[::-1] * y_center[-fade_len:] +
        (1 - weights[::-1]) * right[:fade_len]
    )

    return np.concatenate([left[:-fade_len], y_center, right[fade_len:]])

with rasterio.open(input_dem) as src:
    dem = src.read(1).astype(float)
    dem[dem == src.nodata] = np.nan


dem = dem[:, ::10]
dem = gaussian_filter(dem, sigma=1)
dem = dem - np.nanmin(dem)




row_indices = list(range(0, dem.shape[0], skip))[::-1]
fig, ax = plt.subplots(figsize=(16, 9), dpi=300)

for i, j in enumerate(row_indices):
    row = dem[j]
    if np.all(np.isnan(row)):
        continue

    elevation = np.nan_to_num(row) * vertical_exaggeration

    padded = taper_from_edge_value(
        elevation,
        pad_width=pad_width,
        baseline=baseline,
        edge_jitter=edge_jitter,
        body_jitter=body_jitter,
        fade_len=fade_len
    )

    y = padded + i * line_spacing
    x = np.arange(len(padded))
    ax.plot(x, y, color='white', linewidth=2)

# Joy Division styling
ax.axis('off')
ax.set_facecolor('black')
fig.patch.set_facecolor('black')

# Export
plt.tight_layout()
plt.savefig("joy_plot.png", dpi=300, bbox_inches='tight', pad_inches=0)
# plt.savefig("joy_division_smooth_transition.svg", bbox_inches='tight', pad_inches=0)
plt.show()