Seismic Intensity Map

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LinearSegmentedColormap

# Create grid
x = np.linspace(0, 200, 250)
y = np.linspace(0, 200, 250)
X, Y = np.meshgrid(x, y)

# Epicenter
epi_x, epi_y = 60, 80
R = np.sqrt((X - epi_x)**2 + (Y - epi_y)**2)
R = np.maximum(R, 1)

# Seismic intensity with directional effects
angle = np.arctan2(Y - epi_y, X - epi_x)
directivity = 1 + 0.3 * np.cos(2 * (angle - np.pi/4))
I = 10 * directivity / (1 + R/20) * np.exp(-R/80)
I = np.clip(I, 0, 10)

# Seismic fire theme
plt.style.use('dark_background')
fig, ax = plt.subplots(figsize=(10, 8))
fig.patch.set_facecolor('#100808')
ax.set_facecolor('#100808')

# Custom colormap - seismic fire
colors = ['#100808', '#200a0a', '#400a0a', '#601010', '#802020', 
          '#a04020', '#c06020', '#e08020', '#ffa040', '#ffcc60', '#ffff90']
cmap = LinearSegmentedColormap.from_list('seismic_fire', colors, N=256)

# Filled contour
cf = ax.contourf(X, Y, I, levels=22, cmap=cmap)

# Isoseismal lines
cs = ax.contour(X, Y, I, levels=[2, 4, 6, 8], colors='white', linewidths=0.8, alpha=0.5)
ax.clabel(cs, inline=True, fontsize=9, fmt='MMI %.0f', colors='white')

# Epicenter with glow
ax.scatter([epi_x], [epi_y], color='#ff0000', s=400, marker='*', zorder=6,
           edgecolors='white', linewidths=2)
ax.scatter([epi_x], [epi_y], color='#ff0000', s=800, marker='*', zorder=5, alpha=0.3)

# Fault line
ax.plot([40, 80], [100, 60], '--', color='#ff6666', linewidth=2.5, alpha=0.8)

# Styled colorbar
cbar = plt.colorbar(cf, ax=ax, pad=0.02, shrink=0.85)
cbar.ax.set_facecolor('#100808')
cbar.set_label('Modified Mercalli Intensity', fontsize=11, color='#ffc0a0', labelpad=10)
cbar.ax.yaxis.set_tick_params(color='#ffc0a0')
cbar.outline.set_edgecolor('#402020')
plt.setp(plt.getp(cbar.ax.axes, 'yticklabels'), color='#ffc0a0', fontsize=9)

# Labels
ax.set_xlabel('Distance East (km)', fontsize=13, color='#ffc0a0', fontweight='600', labelpad=10)
ax.set_ylabel('Distance North (km)', fontsize=13, color='#ffc0a0', fontweight='600', labelpad=10)
ax.set_title('Seismic Intensity Map - M6.5 Event', fontsize=16, color='white', 
             fontweight='bold', pad=20)

# Style axes
ax.tick_params(colors='#d0a080', labelsize=10, length=0)
for spine in ax.spines.values():
    spine.set_visible(False)

plt.tight_layout()
plt.show()