{"spec_id":"stereonet-equal-area","library":"letsplot","language":"python","code":"\"\"\" anyplot.ai\nstereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\nLibrary: letsplot 4.10.1 | Python 3.13.13\nQuality: 88/100 | Updated: 2026-06-16\n\"\"\"\n\nimport os\n\nimport numpy as np\nimport pandas as pd\nfrom lets_plot import (\n    LetsPlot,\n    aes,\n    coord_fixed,\n    element_blank,\n    element_rect,\n    element_text,\n    geom_density2d,\n    geom_path,\n    geom_point,\n    geom_polygon,\n    geom_segment,\n    geom_text,\n    ggplot,\n    ggsize,\n    labs,\n    layer_tooltips,\n    scale_color_manual,\n    scale_fill_manual,\n    scale_x_continuous,\n    scale_y_continuous,\n    theme,\n)\nfrom lets_plot.export import ggsave\n\n\nLetsPlot.setup_html()\n\n# Theme-adaptive chrome (Imprint design tokens)\nTHEME = os.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG = \"#FAF8F1\" if THEME == \"light\" else \"#1A1A17\"\nELEVATED_BG = \"#FFFDF6\" if THEME == \"light\" else \"#242420\"\nINK = \"#1A1A17\" if THEME == \"light\" else \"#F0EFE8\"\nINK_SOFT = \"#4A4A44\" if THEME == \"light\" else \"#B8B7B0\"\nINK_MUTED = \"#6B6A63\" if THEME == \"light\" else \"#A8A79F\"\nRULE = \"#D8D4C8\" if THEME == \"light\" else \"#33332E\"\n\n# Imprint categorical palette — first series ALWAYS #009E73\n# 4 abstract structural-feature classes → canonical order 1→4\nIMPRINT_PALETTE = [\"#009E73\", \"#C475FD\", \"#4467A3\", \"#BD8233\"]\n\n# Data - Field measurements from a structural geology mapping campaign\nnp.random.seed(42)\n\n# Bedding planes (NE strike ~045, moderate dip ~35° SE)\nn_bedding = 25\nbedding_strike = np.random.normal(45, 12, n_bedding) % 360\nbedding_dip = np.clip(np.random.normal(35, 8, n_bedding), 5, 85)\n\n# Joint set (N-S striking ~000, steep dip ~78° W)\nn_joints = 20\njoints_strike = np.random.normal(355, 10, n_joints) % 360\njoints_dip = np.clip(np.random.normal(78, 7, n_joints), 10, 89)\n\n# Faults (NW-SE striking ~150, moderate-steep dip ~60°)\nn_faults = 13\nfaults_strike = np.random.uniform(130, 170, n_faults)\nfaults_dip = np.clip(np.random.normal(60, 12, n_faults), 10, 89)\n\n# Foliation (E-W strike ~270, gentle dip ~25° N)\nn_foliation = 12\nfoliation_strike = np.random.normal(270, 8, n_foliation) % 360\nfoliation_dip = np.clip(np.random.normal(25, 6, n_foliation), 5, 60)\n\nstrikes = np.concatenate([bedding_strike, joints_strike, faults_strike, foliation_strike])\ndips = np.concatenate([bedding_dip, joints_dip, faults_dip, foliation_dip])\nfeature_types = [\"Bedding\"] * n_bedding + [\"Joint\"] * n_joints + [\"Fault\"] * n_faults + [\"Foliation\"] * n_foliation\n\n# Equal-area (Schmidt net) lower-hemisphere projection\nstrike_rad = np.radians(strikes)\ndip_rad = np.radians(dips)\ndip_dir_rad = strike_rad + np.pi / 2\n\n# Pole 3D coordinates (lower hemisphere normal to each plane)\npole_x_3d = np.sin(dip_rad) * np.sin(dip_dir_rad)\npole_y_3d = np.sin(dip_rad) * np.cos(dip_dir_rad)\npole_z_3d = -np.cos(dip_rad)\n\n# Lambert azimuthal equal-area projection\npole_scale = 1.0 / np.sqrt(1.0 - pole_z_3d)\npole_x = pole_x_3d * pole_scale\npole_y = pole_y_3d * pole_scale\n\npoles_df = pd.DataFrame(\n    {\n        \"x\": pole_x,\n        \"y\": pole_y,\n        \"feature_type\": feature_types,\n        \"strike\": [f\"{s:.0f}°\" for s in strikes],\n        \"dip\": [f\"{d:.0f}°\" for d in dips],\n    }\n)\n\n# Great circle paths for each plane\ngc_rows = []\nalphas = np.linspace(0, np.pi, 60)\ncos_a = np.cos(alphas)\nsin_a = np.sin(alphas)\n\nfor i in range(len(strikes)):\n    s = strike_rad[i]\n    d = dip_rad[i]\n    dd = dip_dir_rad[i]\n    sv_x, sv_y = np.sin(s), np.cos(s)\n    dv_x = np.cos(d) * np.sin(dd)\n    dv_y = np.cos(d) * np.cos(dd)\n    dv_z = -np.sin(d)\n    px = cos_a * sv_x + sin_a * dv_x\n    py = cos_a * sv_y + sin_a * dv_y\n    pz = sin_a * dv_z\n    sc = 1.0 / np.sqrt(1.0 - pz)\n    gx = px * sc\n    gy = py * sc\n    # Clip to unit circle boundary\n    r2 = gx**2 + gy**2\n    mask_inside = r2 <= 1.0\n    gx_clip = gx[mask_inside]\n    gy_clip = gy[mask_inside]\n    for j in range(len(gx_clip)):\n        gc_rows.append({\"x\": gx_clip[j], \"y\": gy_clip[j], \"group\": i, \"feature_type\": feature_types[i]})\n\ngc_df = pd.DataFrame(gc_rows)\n\n# Primitive circle (outer boundary)\ntheta = np.linspace(0, 2 * np.pi, 200)\ncircle_df = pd.DataFrame({\"x\": np.cos(theta), \"y\": np.sin(theta)})\n\n# Annular mask — covers any density bleed outside the primitive circle with the\n# page background so the stereonet stays cleanly bounded (built as fan quads\n# between r=1 and r=1.6, one filled polygon per angular segment).\nmask_rows = []\nn_seg = 180\nfor k in range(n_seg):\n    a0 = 2 * np.pi * k / n_seg\n    a1 = 2 * np.pi * (k + 1) / n_seg\n    inner0 = (np.cos(a0), np.sin(a0))\n    inner1 = (np.cos(a1), np.sin(a1))\n    outer0 = (1.6 * np.cos(a0), 1.6 * np.sin(a0))\n    outer1 = (1.6 * np.cos(a1), 1.6 * np.sin(a1))\n    for vx, vy in [inner0, outer0, outer1, inner1]:\n        mask_rows.append({\"x\": vx, \"y\": vy, \"group\": k})\nmask_df = pd.DataFrame(mask_rows)\n\n# Tick marks every 10° (longer every 30°)\ntick_rows = []\nfor deg in range(0, 360, 10):\n    rad = np.radians(deg)\n    inner_r = 0.95 if deg % 30 != 0 else 0.92\n    tick_rows.append({\"x\": np.sin(rad) * inner_r, \"y\": np.cos(rad) * inner_r, \"xend\": np.sin(rad), \"yend\": np.cos(rad)})\ntick_df = pd.DataFrame(tick_rows)\n\n# Reference cross lines (N-S, E-W)\nref_df = pd.DataFrame(\n    [{\"x\": 0, \"y\": -1, \"xend\": 0, \"yend\": 1, \"group\": 0}, {\"x\": -1, \"y\": 0, \"xend\": 1, \"yend\": 0, \"group\": 1}]\n)\n\n# Cardinal direction labels\ncardinal_positions = [(0, 1.14, \"N\"), (1.14, 0, \"E\"), (0, -1.14, \"S\"), (-1.14, 0, \"W\")]\nlabel_df = pd.DataFrame(\n    {\n        \"x\": [c[0] for c in cardinal_positions],\n        \"y\": [c[1] for c in cardinal_positions],\n        \"label\": [c[2] for c in cardinal_positions],\n    }\n)\n\n# Mean pole annotations and confidence ellipses per feature type\nmean_annotations = []\nellipse_rows = []\nfor idx, ft in enumerate([\"Bedding\", \"Joint\", \"Fault\", \"Foliation\"]):\n    mask = np.array(feature_types) == ft\n    mx = pole_x[mask].mean()\n    my = pole_y[mask].mean()\n    # Circular mean for strike (handles wraparound at 0/360°)\n    s_rad = np.radians(strikes[mask])\n    ms = np.degrees(np.arctan2(np.sin(s_rad).mean(), np.cos(s_rad).mean())) % 360\n    md = dips[mask].mean()\n    # Smart label nudge: push outward from cluster centroid, away from cardinal labels\n    nudge_x, nudge_y = 0.0, -0.20\n    for cx, cy, _ in cardinal_positions:\n        dist = np.sqrt((mx - cx) ** 2 + (my - cy) ** 2)\n        if dist < 0.45:\n            # Push label away from the nearby cardinal direction, toward the centre\n            dx, dy = mx - cx, my - cy\n            norm = max(np.sqrt(dx**2 + dy**2), 0.01)\n            nudge_x = dx / norm * 0.20\n            nudge_y = dy / norm * 0.20 - 0.06\n    mean_annotations.append(\n        {\"x\": mx, \"y\": my, \"label\": f\"{ft}\\n{ms:.0f}/{md:.0f}\", \"nudge_x\": nudge_x, \"nudge_y\": nudge_y}\n    )\n    # Confidence ellipse (1-sigma) around each cluster\n    px_ft = pole_x[mask]\n    py_ft = pole_y[mask]\n    cov = np.cov(px_ft, py_ft)\n    eigvals, eigvecs = np.linalg.eigh(cov)\n    angle = np.arctan2(eigvecs[1, 1], eigvecs[0, 1])\n    t = np.linspace(0, 2 * np.pi, 40)\n    ex = np.sqrt(eigvals[1]) * np.cos(t)\n    ey = np.sqrt(eigvals[0]) * np.sin(t)\n    rx = mx + ex * np.cos(angle) - ey * np.sin(angle)\n    ry = my + ex * np.sin(angle) + ey * np.cos(angle)\n    for j in range(len(t)):\n        ellipse_rows.append({\"x\": rx[j], \"y\": ry[j], \"group\": idx, \"feature_type\": ft})\n\nellipse_df = pd.DataFrame(ellipse_rows)\nmean_df = pd.DataFrame(mean_annotations)\n\n# Tooltips for poles showing strike/dip\npole_tooltips = layer_tooltips().line(\"@feature_type\").line(\"Strike: @strike\").line(\"Dip: @dip\")\n\n# Build individual mean label layers to handle per-label nudge offsets\nmean_label_layers = []\nfor _, row in mean_df.iterrows():\n    single_df = pd.DataFrame([{\"x\": row[\"x\"], \"y\": row[\"y\"], \"label\": row[\"label\"]}])\n    mean_label_layers.append(\n        geom_text(\n            aes(x=\"x\", y=\"y\", label=\"label\"),\n            data=single_df,\n            size=7,\n            color=INK,\n            nudge_x=row[\"nudge_x\"],\n            nudge_y=row[\"nudge_y\"],\n            fontface=\"italic\",\n            show_legend=False,\n        )\n    )\n\n# Plot\nplot = (\n    ggplot()\n    # Reference cross lines\n    + geom_segment(aes(x=\"x\", y=\"y\", xend=\"xend\", yend=\"yend\"), data=ref_df, color=RULE, size=0.5, linetype=\"dashed\")\n    # Density contours (Kamb-style preferred-orientation clustering)\n    + geom_density2d(\n        aes(x=\"x\", y=\"y\"),\n        data=poles_df,\n        color=INK_MUTED,\n        alpha=0.55,\n        size=0.6,\n        bins=6,\n        kernel=\"gaussian\",\n        adjust=0.8,\n        show_legend=False,\n    )\n    # Great circles\n    + geom_path(aes(x=\"x\", y=\"y\", group=\"group\", color=\"feature_type\"), data=gc_df, size=0.35, alpha=0.18)\n    # Confidence ellipses around each cluster\n    + geom_polygon(\n        aes(x=\"x\", y=\"y\", group=\"group\", fill=\"feature_type\"), data=ellipse_df, alpha=0.12, size=0, show_legend=False\n    )\n    # Poles to planes with tooltips\n    + geom_point(aes(x=\"x\", y=\"y\", color=\"feature_type\"), data=poles_df, size=4.5, alpha=0.88, tooltips=pole_tooltips)\n    # Mean pole markers (diamond shape)\n    + geom_point(aes(x=\"x\", y=\"y\"), data=mean_df, size=9, shape=18, color=INK, alpha=0.95, show_legend=False)\n)\n\n# Add per-label mean annotation layers\nfor layer in mean_label_layers:\n    plot = plot + layer\n\nplot = (\n    plot\n    # Annular mask hides density bleed outside the primitive circle\n    + geom_polygon(aes(x=\"x\", y=\"y\", group=\"group\"), data=mask_df, fill=PAGE_BG, color=PAGE_BG, size=0)\n    # Primitive circle\n    + geom_path(aes(x=\"x\", y=\"y\"), data=circle_df, color=INK, size=1.4)\n    # Tick marks\n    + geom_segment(aes(x=\"x\", y=\"y\", xend=\"xend\", yend=\"yend\"), data=tick_df, color=INK, size=0.9)\n    # Cardinal labels\n    + geom_text(aes(x=\"x\", y=\"y\", label=\"label\"), data=label_df, size=9, color=INK, fontface=\"bold\")\n    # Color scale (Imprint palette)\n    + scale_color_manual(values=IMPRINT_PALETTE, name=\"Feature Type\")\n    + scale_fill_manual(values=IMPRINT_PALETTE)\n    + coord_fixed()\n    + scale_x_continuous(limits=[-1.32, 1.32])\n    + scale_y_continuous(limits=[-1.32, 1.32])\n    + labs(title=\"stereonet-equal-area · python · letsplot · anyplot.ai\")\n    + theme(\n        plot_title=element_text(size=18, hjust=0.5, face=\"bold\", color=INK, margin=[0, 0, 12, 0]),\n        legend_title=element_text(size=14, face=\"bold\", color=INK),\n        legend_text=element_text(size=12, color=INK_SOFT),\n        legend_position=[0.86, 0.16],\n        legend_justification=[0.5, 0.5],\n        legend_background=element_rect(fill=ELEVATED_BG, color=INK_SOFT, size=0.5),\n        axis_title=element_blank(),\n        axis_text=element_blank(),\n        axis_ticks=element_blank(),\n        axis_line=element_blank(),\n        panel_grid=element_blank(),\n        plot_background=element_rect(fill=PAGE_BG, color=PAGE_BG),\n        panel_background=element_rect(fill=PAGE_BG, color=PAGE_BG),\n    )\n    + ggsize(600, 600)\n)\n\nggsave(plot, f\"plot-{THEME}.png\", path=\".\", scale=4)\nggsave(plot, f\"plot-{THEME}.html\", path=\".\")\n"}