{"spec_id":"stereonet-equal-area","library":"plotly","language":"python","code":"\"\"\" anyplot.ai\nstereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\nLibrary: plotly 6.8.0 | Python 3.13.13\nQuality: 92/100 | Updated: 2026-06-16\n\"\"\"\n\nimport os\n\nimport numpy as np\nimport plotly.graph_objects as go\nfrom scipy.stats import gaussian_kde\n\n\n# Theme-adaptive chrome (Imprint palette)\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 = \"rgba(26,26,23,0.15)\" if THEME == \"light\" else \"rgba(240,239,232,0.15)\"\n\n# Imprint categorical palette — feature types are abstract, so canonical order.\n# Bedding takes the brand green (always-first series).\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-striking, moderate SE dip\nbedding_strike = np.random.normal(45, 8, 20)\nbedding_dip = np.random.normal(30, 5, 20)\n\n# Joint set 1: N-S striking, steep E dip\njoint1_strike = np.random.normal(0, 10, 15)\njoint1_dip = np.random.normal(80, 5, 15)\n\n# Joint set 2: E-W striking, steep S dip\njoint2_strike = np.random.normal(90, 12, 12)\njoint2_dip = np.random.normal(75, 8, 12)\n\n# Faults: NW-SE striking, moderate dip\nfault_strike = np.random.normal(130, 15, 8)\nfault_dip = np.random.normal(55, 10, 8)\n\nstrikes = np.concatenate([bedding_strike, joint1_strike, joint2_strike, fault_strike])\ndips = np.concatenate([bedding_dip, joint1_dip, joint2_dip, fault_dip])\nfeature_types = [\"Bedding\"] * 20 + [\"Joint Set 1\"] * 15 + [\"Joint Set 2\"] * 12 + [\"Fault\"] * 8\nstrikes = strikes % 360\ndips = np.clip(dips, 0, 90)\n\ntype_colors = {\n    \"Bedding\": IMPRINT_PALETTE[0],\n    \"Joint Set 1\": IMPRINT_PALETTE[1],\n    \"Joint Set 2\": IMPRINT_PALETTE[2],\n    \"Fault\": IMPRINT_PALETTE[3],\n}\n\n# Equal-area (Schmidt) projection of poles to planes\n# Pole to plane: trend = dip direction = strike + 90°, plunge = 90° - dip\npole_trend_rad = np.radians((strikes + 90) % 360)\npole_plunge_rad = np.radians(90 - dips)\npole_r = np.sqrt(2) * np.sin((np.pi / 2 - pole_plunge_rad) / 2)\npole_x = pole_r * np.sin(pole_trend_rad)\npole_y = pole_r * np.cos(pole_trend_rad)\n\n# Plot\nfig = go.Figure()\n\n# Density contours on pole data (Kamb-style KDE).\n# Continuous density → single-polarity imprint_seq (green → blue), with rising alpha\n# so low density fades into the page and the peak reads clearly over the markers.\nxy_poles = np.vstack([pole_x, pole_y])\nkde = gaussian_kde(xy_poles, bw_method=0.2)\ngrid_n = 150\ngx = np.linspace(-1.02, 1.02, grid_n)\ngy = np.linspace(-1.02, 1.02, grid_n)\nGX, GY = np.meshgrid(gx, gy)\nZ = kde(np.vstack([GX.ravel(), GY.ravel()])).reshape(grid_n, grid_n)\nZ[GX**2 + GY**2 > 1.0] = np.nan\n\ndensity_colorscale = [\n    [0.0, \"rgba(0,158,115,0.0)\"],  # imprint_seq low — transparent green\n    [0.25, \"rgba(0,158,115,0.16)\"],\n    [0.5, \"rgba(0,158,115,0.34)\"],\n    [0.75, \"rgba(68,103,163,0.50)\"],\n    [1.0, \"rgba(68,103,163,0.66)\"],  # imprint_seq high — blue\n]\n\nfig.add_trace(\n    go.Contour(\n        x=gx,\n        y=gy,\n        z=Z,\n        colorscale=density_colorscale,\n        showscale=False,\n        contours={\"coloring\": \"fill\", \"showlines\": True, \"showlabels\": False},\n        line={\"color\": RULE, \"width\": 0.8},\n        ncontours=6,\n        showlegend=False,\n        hoverinfo=\"skip\",\n    )\n)\n\n# Stereonet grid - primitive circle\ntheta_circ = np.linspace(0, 2 * np.pi, 361)\nfig.add_trace(\n    go.Scatter(\n        x=np.cos(theta_circ).tolist(),\n        y=np.sin(theta_circ).tolist(),\n        mode=\"lines\",\n        line={\"color\": INK, \"width\": 2.2},\n        showlegend=False,\n        hoverinfo=\"skip\",\n    )\n)\n\n# Grid: small circles at 30° and 60° inclination from center\ngrid_x, grid_y = [], []\nfor inc_deg in [30, 60]:\n    r_grid = np.sqrt(2) * np.sin(np.radians(inc_deg) / 2)\n    grid_x.extend(list(r_grid * np.cos(theta_circ)) + [None])\n    grid_y.extend(list(r_grid * np.sin(theta_circ)) + [None])\n\n# Grid: N-S and E-W diameter lines\nfor angle in [0, np.pi / 2]:\n    grid_x.extend([-np.sin(angle), np.sin(angle), None])\n    grid_y.extend([-np.cos(angle), np.cos(angle), None])\n\nfig.add_trace(\n    go.Scatter(\n        x=grid_x,\n        y=grid_y,\n        mode=\"lines\",\n        line={\"color\": RULE, \"width\": 0.8, \"dash\": \"dot\"},\n        showlegend=False,\n        hoverinfo=\"skip\",\n    )\n)\n\n# Tick marks every 10° around perimeter\ntick_x, tick_y = [], []\nfor tick_deg in range(0, 360, 10):\n    tick_rad = np.radians(tick_deg)\n    tick_len = 0.95 if tick_deg % 30 != 0 else 0.93\n    tick_x.extend([np.sin(tick_rad), tick_len * np.sin(tick_rad), None])\n    tick_y.extend([np.cos(tick_rad), tick_len * np.cos(tick_rad), None])\n\nfig.add_trace(\n    go.Scatter(\n        x=tick_x, y=tick_y, mode=\"lines\", line={\"color\": INK_SOFT, \"width\": 1.3}, showlegend=False, hoverinfo=\"skip\"\n    )\n)\n\n# Great circles for representative planes (subset to avoid clutter)\ngc_alpha = np.linspace(0, np.pi, 180)\ngc_indices = list(range(0, 20, 5)) + list(range(20, 35, 7)) + list(range(35, 47, 6)) + list(range(47, 55, 4))\n\nfor idx in gc_indices:\n    trend_rad = np.radians((strikes[idx] + 90) % 360)\n    plunge_rad = np.radians(90 - dips[idx])\n    cos_p, sin_p = np.cos(plunge_rad), np.sin(plunge_rad)\n    cos_t, sin_t = np.cos(trend_rad), np.sin(trend_rad)\n\n    # v1: horizontal vector perpendicular to pole trend\n    v1 = np.array([-cos_t, sin_t, 0.0])\n    # v2 = pole_vector × v1\n    v2 = np.array([-sin_p * sin_t, -sin_p * cos_t, cos_p])\n\n    # Great circle parameterization (lower hemisphere: alpha in [0, pi])\n    gc_pts = np.outer(np.cos(gc_alpha), v1) + np.outer(np.sin(gc_alpha), v2)\n    gc_vz = gc_pts[:, 2]\n    lower = gc_vz >= 0\n    if lower.sum() < 2:\n        continue\n\n    gc_vx, gc_vy, gc_vz = gc_pts[lower, 0], gc_pts[lower, 1], gc_pts[lower, 2]\n    gc_plunge = np.arcsin(np.clip(gc_vz, -1, 1))\n    gc_trend = np.arctan2(gc_vx, gc_vy)\n    gc_r = np.sqrt(2) * np.sin((np.pi / 2 - gc_plunge) / 2)\n    gc_proj_x = gc_r * np.sin(gc_trend)\n    gc_proj_y = gc_r * np.cos(gc_trend)\n\n    fig.add_trace(\n        go.Scatter(\n            x=gc_proj_x.tolist(),\n            y=gc_proj_y.tolist(),\n            mode=\"lines\",\n            line={\"color\": type_colors[feature_types[idx]], \"width\": 2.0},\n            opacity=0.65,\n            legendgroup=feature_types[idx],\n            showlegend=False,\n            hoverinfo=\"skip\",\n        )\n    )\n\n# Plot poles by feature type with Plotly customdata + hovertemplate.\n# Marker edge matches the page bg for a clean halo that separates overlapping poles.\nfor feat_type in [\"Bedding\", \"Joint Set 1\", \"Joint Set 2\", \"Fault\"]:\n    mask = np.array([t == feat_type for t in feature_types])\n    customdata = np.column_stack([strikes[mask], dips[mask], (strikes[mask] + 90) % 360])\n    fig.add_trace(\n        go.Scatter(\n            x=pole_x[mask].tolist(),\n            y=pole_y[mask].tolist(),\n            mode=\"markers\",\n            name=feat_type,\n            legendgroup=feat_type,\n            marker={\n                \"size\": 12,\n                \"color\": type_colors[feat_type],\n                \"line\": {\"width\": 1.5, \"color\": PAGE_BG},\n                \"symbol\": \"circle\",\n            },\n            customdata=customdata,\n            hovertemplate=(\n                f\"<b>{feat_type}</b><br>\"\n                \"Strike: %{customdata[0]:.0f}°<br>\"\n                \"Dip: %{customdata[1]:.0f}°<br>\"\n                \"Dip Direction: %{customdata[2]:.0f}°\"\n                \"<extra></extra>\"\n            ),\n        )\n    )\n\n# Cardinal direction and degree labels\nfor label, lx, ly in [(\"N\", 0, 1.12), (\"E\", 1.12, 0), (\"S\", 0, -1.12), (\"W\", -1.12, 0)]:\n    fig.add_annotation(\n        x=lx,\n        y=ly,\n        text=f\"<b>{label}</b>\",\n        showarrow=False,\n        font={\"size\": 17, \"color\": INK},\n        xanchor=\"center\",\n        yanchor=\"middle\",\n    )\n\nfor deg in range(0, 360, 30):\n    if deg % 90 == 0:\n        continue\n    rad = np.radians(deg)\n    fig.add_annotation(\n        x=1.075 * np.sin(rad),\n        y=1.075 * np.cos(rad),\n        text=f\"{deg}°\",\n        showarrow=False,\n        font={\"size\": 10, \"color\": INK_MUTED},\n        xanchor=\"center\",\n        yanchor=\"middle\",\n    )\n\n# Interactive buttons for density contour toggle (Plotly-specific feature)\nn_traces = len(fig.data)\ndensity_visible_on = [True] * n_traces\ndensity_visible_off = [True] * n_traces\ndensity_visible_off[0] = False  # First trace is the density contour\n\n# Style\nfig.update_layout(\n    autosize=False,\n    width=600,\n    height=600,\n    updatemenus=[\n        {\n            \"type\": \"buttons\",\n            \"direction\": \"left\",\n            \"buttons\": [\n                {\"label\": \"Show Density\", \"method\": \"update\", \"args\": [{\"visible\": density_visible_on}]},\n                {\"label\": \"Hide Density\", \"method\": \"update\", \"args\": [{\"visible\": density_visible_off}]},\n            ],\n            \"x\": 0.01,\n            \"y\": 0.04,\n            \"xanchor\": \"left\",\n            \"yanchor\": \"bottom\",\n            \"bgcolor\": ELEVATED_BG,\n            \"bordercolor\": INK_SOFT,\n            \"font\": {\"size\": 11, \"color\": INK},\n        }\n    ],\n    title={\n        \"text\": (\n            \"stereonet-equal-area · python · plotly · anyplot.ai\"\n            f\"<br><sup style='color:{INK_SOFT}'>Lower Hemisphere, Equal-Area (Schmidt) Projection</sup>\"\n        ),\n        \"font\": {\"size\": 18, \"color\": INK},\n        \"x\": 0.5,\n        \"xanchor\": \"center\",\n        \"y\": 0.97,\n        \"yanchor\": \"top\",\n    },\n    xaxis={\n        \"scaleanchor\": \"y\",\n        \"scaleratio\": 1,\n        \"showgrid\": False,\n        \"zeroline\": False,\n        \"showticklabels\": False,\n        \"showline\": False,\n        \"range\": [-1.2, 1.2],\n    },\n    yaxis={\"showgrid\": False, \"zeroline\": False, \"showticklabels\": False, \"showline\": False, \"range\": [-1.2, 1.2]},\n    legend={\n        \"title\": {\"text\": \"Feature Type\", \"font\": {\"size\": 13, \"color\": INK}},\n        \"font\": {\"size\": 12, \"color\": INK_SOFT},\n        \"bgcolor\": ELEVATED_BG,\n        \"bordercolor\": INK_SOFT,\n        \"borderwidth\": 1,\n        \"x\": 0.99,\n        \"y\": 0.99,\n        \"xanchor\": \"right\",\n        \"yanchor\": \"top\",\n    },\n    paper_bgcolor=PAGE_BG,\n    plot_bgcolor=PAGE_BG,\n    font={\"color\": INK},\n    margin={\"l\": 20, \"r\": 20, \"t\": 70, \"b\": 20},\n)\n\n# Annotation highlighting dominant bedding cluster\nbedding_pole_x = pole_x[:20].mean()\nbedding_pole_y = pole_y[:20].mean()\nfig.add_annotation(\n    x=bedding_pole_x + 0.28,\n    y=bedding_pole_y - 0.18,\n    ax=bedding_pole_x,\n    ay=bedding_pole_y,\n    text=\"Dominant NE-striking<br>bedding fabric\",\n    showarrow=True,\n    arrowhead=2,\n    arrowsize=1,\n    arrowwidth=1.5,\n    arrowcolor=IMPRINT_PALETTE[0],\n    font={\"size\": 11, \"color\": INK},\n    bgcolor=ELEVATED_BG,\n    bordercolor=IMPRINT_PALETTE[0],\n    borderwidth=1,\n    borderpad=4,\n)\n\n# Save — square canvas (radially symmetric stereonet): 600×600 × scale 4 = 2400×2400\nfig.write_image(f\"plot-{THEME}.png\", width=600, height=600, scale=4)\nfig.write_html(f\"plot-{THEME}.html\", include_plotlyjs=\"cdn\")\n"}