{"spec_id":"stereonet-equal-area","library":"altair","language":"python","code":"\"\"\" anyplot.ai\nstereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\nLibrary: altair 6.2.1 | Python 3.13.13\nQuality: 90/100 | Updated: 2026-06-16\n\"\"\"\n\nimport os\n\nimport altair as alt\nimport numpy as np\nimport pandas as pd\nfrom contourpy import contour_generator\n\n\n# Theme-adaptive chrome\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\"\n\n# Imprint palette — categorical, theme-independent, hybrid-v3 sort order.\n# Geological feature types are abstract categories → canonical order 1→3.\nIMPRINT = [\"#009E73\", \"#C475FD\", \"#4467A3\"]\nFEATURES = [\"Bedding\", \"Fault\", \"Joint\"]\ncolor_map = dict(zip(FEATURES, IMPRINT, strict=True))\ncolor_scale = alt.Scale(domain=FEATURES, range=IMPRINT)\n# Continuous density → imprint_seq (single-polarity, brand green → blue).\nIMPRINT_SEQ = [\"#009E73\", \"#4467A3\"]\n\n# Data — field measurements of geological structures\nnp.random.seed(42)\n\nbedding_strike = np.random.normal(45, 12, 25)\nbedding_dip = np.random.normal(35, 8, 25)\n\nfault_strike = np.random.normal(310, 15, 18)\nfault_dip = np.random.normal(70, 10, 18)\n\njoint_strike = np.random.normal(90, 10, 22)\njoint_dip = np.random.normal(80, 7, 22)\n\nstrikes = np.concatenate([bedding_strike, fault_strike, joint_strike]) % 360\ndips = np.clip(np.concatenate([bedding_dip, fault_dip, joint_dip]), 0, 90)\nfeature_types = [\"Bedding\"] * len(bedding_strike) + [\"Fault\"] * len(fault_strike) + [\"Joint\"] * len(joint_strike)\n\n# Equal-area (Schmidt) projection: poles to planes\npole_trend = np.radians((strikes + 90) % 360)\npole_r = np.sqrt(2) * np.sin(np.radians(dips) / 2)\npole_x = pole_r * np.sin(pole_trend)\npole_y = pole_r * np.cos(pole_trend)\n\ndf_poles = pd.DataFrame(\n    {\"x\": pole_x, \"y\": pole_y, \"feature_type\": feature_types, \"strike\": np.round(strikes, 1), \"dip\": np.round(dips, 1)}\n)\n\n# Great circles for each measurement\ngc_rows = []\nfor i in range(len(strikes)):\n    s_rad = np.radians(strikes[i])\n    d_rad = np.radians(dips[i])\n    dd_rad = s_rad + np.pi / 2\n    v1 = np.array([np.sin(s_rad), np.cos(s_rad), 0.0])\n    v2 = np.array([np.cos(d_rad) * np.sin(dd_rad), np.cos(d_rad) * np.cos(dd_rad), -np.sin(d_rad)])\n    for j, rake in enumerate(np.linspace(0, np.pi, 61)):\n        line = np.cos(rake) * v1 + np.sin(rake) * v2\n        if line[2] > 0:\n            line = -line\n        plunge = np.arcsin(-line[2])\n        trend = np.arctan2(line[0], line[1])\n        r = np.sqrt(2) * np.sin((np.pi / 2 - plunge) / 2)\n        gc_rows.append(\n            {\"x\": r * np.sin(trend), \"y\": r * np.cos(trend), \"feature_type\": feature_types[i], \"gc_id\": i, \"order\": j}\n        )\ndf_gc = pd.DataFrame(gc_rows)\n\n# Primitive circle\nr_prim = np.sqrt(2) * np.sin(np.pi / 4)\ntheta_circ = np.linspace(0, 2 * np.pi, 361)\ndf_circle = pd.DataFrame({\"x\": r_prim * np.sin(theta_circ), \"y\": r_prim * np.cos(theta_circ), \"order\": range(361)})\n\n# Tick marks every 10 degrees around the perimeter (longer at 30-degree marks)\ntick_rows = []\nfor deg in range(0, 360, 10):\n    rad = np.radians(deg)\n    tick_len = 0.06 if deg % 30 == 0 else 0.04\n    tick_rows.append({\"x\": r_prim * np.sin(rad), \"y\": r_prim * np.cos(rad), \"tid\": deg, \"order\": 0})\n    tick_rows.append(\n        {\"x\": (r_prim - tick_len) * np.sin(rad), \"y\": (r_prim - tick_len) * np.cos(rad), \"tid\": deg, \"order\": 1}\n    )\ndf_ticks = pd.DataFrame(tick_rows)\n\n# Cardinal labels\nlbl_r = r_prim + 0.09\ndf_dirs = pd.DataFrame(\n    {\"x\": [0, lbl_r, 0, -lbl_r], \"y\": [lbl_r + 0.03, 0, -lbl_r - 0.03, 0], \"label\": [\"N\", \"E\", \"S\", \"W\"]}\n)\n\n# Equal-area net grid circles at 30 and 60 degree dip\ngrid_rows = []\nfor dip_g in [30, 60]:\n    r_g = np.sqrt(2) * np.sin(np.radians(dip_g) / 2)\n    for j, t in enumerate(np.linspace(0, 2 * np.pi, 181)):\n        grid_rows.append({\"x\": r_g * np.sin(t), \"y\": r_g * np.cos(t), \"level\": dip_g, \"order\": j})\ndf_grid = pd.DataFrame(grid_rows)\n\n# Grid cross lines (N-S, E-W)\ndf_cross = pd.DataFrame(\n    {\n        \"x\": [0, 0, -r_prim, r_prim],\n        \"y\": [-r_prim, r_prim, 0, 0],\n        \"line_id\": [\"NS\", \"NS\", \"EW\", \"EW\"],\n        \"order\": [0, 1, 0, 1],\n    }\n)\n\n# Density field via the Kamb (1959) counting-circle statistic — the standard\n# structural-geology method. A counting circle of fixed angular radius is sized\n# so that, for a uniform (random) fabric, the expected count within it has a\n# standard deviation of `sigma`; the local density is then the number of poles\n# falling inside the circle, expressed in units of sigma above that uniform\n# expectation. Contours at integer-sigma levels mark statistically significant\n# preferred orientations.\nn_grid = 160\ngx = np.linspace(-r_prim, r_prim, n_grid)\ngy = np.linspace(-r_prim, r_prim, n_grid)\ngxx, gyy = np.meshgrid(gx, gy)\n\n# Lower-hemisphere 3-D unit vectors for the poles, recovered from their\n# equal-area positions (angle from the vertical = 2*arcsin(r/sqrt(2))).\npole_th = 2.0 * np.arcsin(np.clip(pole_r / np.sqrt(2), 0.0, 1.0))\npole_vec = np.column_stack(\n    [np.sin(pole_th) * np.sin(pole_trend), np.sin(pole_th) * np.cos(pole_trend), -np.cos(pole_th)]\n)\n\n# Grid nodes mapped to the same lower-hemisphere unit vectors.\nnode_r = np.hypot(gxx, gyy)\nnode_th = 2.0 * np.arcsin(np.clip(node_r / np.sqrt(2), 0.0, 1.0))\nnode_tr = np.arctan2(gxx, gyy)\nnode_vec = np.stack([np.sin(node_th) * np.sin(node_tr), np.sin(node_th) * np.cos(node_tr), -np.cos(node_th)], axis=-1)\n\n# Kamb counting circle: cos-distance threshold and the sigma normalisation.\nn_poles = len(pole_vec)\nsigma = 3.0\nkamb_cos = 1.0 - sigma**2 / (n_poles + sigma**2)\nkamb_units = np.sqrt(n_poles * kamb_cos * (1.0 - kamb_cos))\n\n# Count poles inside each node's counting circle (cos of angular distance >=\n# threshold), then convert to sigma units with a 0.5 continuity correction.\ncos_dist = node_vec @ pole_vec.T\ncounts = (cos_dist >= kamb_cos).sum(axis=-1)\ndensity = np.clip((counts - 0.5) / kamb_units, 0.0, None)\ndensity[gxx**2 + gyy**2 > r_prim**2] = 0.0\n\n# Smooth, continuous contour polylines via contourpy (the same iso-path engine\n# matplotlib uses) — connected vertex chains, not the disconnected per-cell\n# segments a marching-squares hand-roll would produce. Contour at integer-sigma\n# levels (2σ, 4σ, …); fall back to relative levels for a near-uniform fabric.\ndmax = float(np.nanmax(density))\ncontour_levels = np.arange(2.0, dmax, 2.0)\nif len(contour_levels) < 2:\n    contour_levels = np.linspace(dmax * 0.3, dmax * 0.9, 4)\ncgen = contour_generator(gx, gy, density, line_type=\"Separate\")\ncontour_rows = []\nfor li, level in enumerate(contour_levels):\n    for si, seg in enumerate(cgen.lines(float(level))):\n        for oi, (cx, cy) in enumerate(seg):\n            if cx * cx + cy * cy <= r_prim**2 + 1e-9:\n                contour_rows.append({\"x\": cx, \"y\": cy, \"seg_id\": f\"{li}_{si}\", \"level\": float(li), \"order\": oi})\ndf_contours = (\n    pd.DataFrame(contour_rows) if contour_rows else pd.DataFrame(columns=[\"x\", \"y\", \"seg_id\", \"level\", \"order\"])\n)\n\n# Mean orientation per feature type; labels pushed radially outward so they\n# clear the dense pole clusters they annotate.\nmean_rows = []\nfor ft in FEATURES:\n    mask = np.array(feature_types) == ft\n    mx, my = float(np.mean(pole_x[mask])), float(np.mean(pole_y[mask]))\n    ms, md = float(np.mean(strikes[mask])), float(np.mean(dips[mask]))\n    # Offset the label radially off the cluster centroid, but cap its radius so\n    # near-perimeter clusters don't push the text onto the cardinal labels.\n    norm = np.hypot(mx, my)\n    if norm > 1e-6:\n        lr = min(norm + 0.30, 0.86)\n        lx, ly = mx / norm * lr, my / norm * lr\n    else:\n        lx, ly = mx, my + 0.30\n    mean_rows.append({\"x\": mx, \"y\": my, \"lx\": lx, \"ly\": ly, \"feature_type\": ft, \"label\": f\"μ {ms:.0f}°/{md:.0f}°\"})\ndf_means = pd.DataFrame(mean_rows)\n\n# Plot\nx_enc = alt.X(\"x:Q\", axis=None, scale=alt.Scale(domain=[-1.3, 1.3]))\ny_enc = alt.Y(\"y:Q\", axis=None, scale=alt.Scale(domain=[-1.3, 1.3]))\n\n# Interactive legend selection for highlighting by feature type\nselection = alt.selection_point(fields=[\"feature_type\"], bind=\"legend\")\n\n# Density contour lines (imprint_seq, drawn beneath the data)\ncontour_layer = (\n    alt.Chart(df_contours)\n    .mark_line(strokeWidth=2.4, strokeCap=\"round\", opacity=0.85)\n    .encode(\n        x=x_enc,\n        y=y_enc,\n        detail=\"seg_id:N\",\n        order=\"order:O\",\n        color=alt.Color(\n            \"level:Q\", scale=alt.Scale(range=IMPRINT_SEQ, domain=[0, len(contour_levels) - 1]), legend=None\n        ),\n    )\n    if len(df_contours) > 0\n    else alt.Chart(pd.DataFrame({\"x\": [0], \"y\": [0]})).mark_point(size=0).encode(x=\"x:Q\", y=\"y:Q\")\n)\n\ngrid_circles = (\n    alt.Chart(df_grid)\n    .mark_line(strokeWidth=0.8, color=INK_SOFT, opacity=0.3)\n    .encode(x=x_enc, y=y_enc, detail=\"level:N\", order=\"order:O\")\n)\n\ncross_lines = (\n    alt.Chart(df_cross)\n    .mark_line(strokeWidth=0.8, color=INK_SOFT, opacity=0.3)\n    .encode(x=x_enc, y=y_enc, detail=\"line_id:N\", order=\"order:O\")\n)\n\ngreat_circles = (\n    alt.Chart(df_gc)\n    .mark_line(strokeWidth=1.1)\n    .encode(\n        x=x_enc,\n        y=y_enc,\n        detail=\"gc_id:N\",\n        order=\"order:O\",\n        color=alt.Color(\"feature_type:N\", scale=color_scale, legend=None),\n        opacity=alt.condition(selection, alt.value(0.32), alt.value(0.05)),\n    )\n    .add_params(selection)\n)\n\nprim_circle = alt.Chart(df_circle).mark_line(strokeWidth=2.5, color=INK).encode(x=x_enc, y=y_enc, order=\"order:O\")\n\ntick_marks = (\n    alt.Chart(df_ticks).mark_line(strokeWidth=1.4, color=INK).encode(x=x_enc, y=y_enc, detail=\"tid:N\", order=\"order:O\")\n)\n\ndir_labels = (\n    alt.Chart(df_dirs).mark_text(fontSize=15, fontWeight=\"bold\", color=INK).encode(x=x_enc, y=y_enc, text=\"label:N\")\n)\n\npoles_layer = (\n    alt.Chart(df_poles)\n    .mark_point(filled=True, size=130, stroke=PAGE_BG, strokeWidth=1.1)\n    .encode(\n        x=x_enc,\n        y=y_enc,\n        color=alt.Color(\n            \"feature_type:N\",\n            scale=color_scale,\n            title=\"Feature Type\",\n            legend=alt.Legend(\n                titleFontSize=12,\n                labelFontSize=11,\n                symbolSize=140,\n                orient=\"top-right\",\n                titleColor=INK,\n                labelColor=INK_SOFT,\n                fillColor=ELEVATED_BG,\n                strokeColor=INK_SOFT,\n                padding=8,\n                cornerRadius=4,\n            ),\n        ),\n        opacity=alt.condition(selection, alt.value(0.92), alt.value(0.12)),\n        tooltip=[\n            alt.Tooltip(\"feature_type:N\", title=\"Type\"),\n            alt.Tooltip(\"strike:Q\", title=\"Strike (°)\"),\n            alt.Tooltip(\"dip:Q\", title=\"Dip (°)\"),\n        ],\n    )\n    .add_params(selection)\n)\n\n# Mean orientation cross markers\nmean_markers = (\n    alt.Chart(df_means)\n    .mark_point(shape=\"cross\", size=260, strokeWidth=3)\n    .encode(\n        x=x_enc,\n        y=y_enc,\n        color=alt.Color(\"feature_type:N\", scale=color_scale, legend=None),\n        opacity=alt.condition(selection, alt.value(1.0), alt.value(0.12)),\n    )\n    .add_params(selection)\n)\n\nmean_labels = (\n    alt.Chart(df_means)\n    .mark_text(fontSize=11, fontWeight=\"bold\", color=INK)\n    .encode(\n        x=alt.X(\"lx:Q\", axis=None, scale=alt.Scale(domain=[-1.3, 1.3])),\n        y=alt.Y(\"ly:Q\", axis=None, scale=alt.Scale(domain=[-1.3, 1.3])),\n        text=\"label:N\",\n        opacity=alt.condition(selection, alt.value(0.95), alt.value(0.1)),\n    )\n    .add_params(selection)\n)\n\nTITLE = \"stereonet-equal-area · python · altair · anyplot.ai\"\ntitle_fontsize = round(16 * 67 / len(TITLE)) if len(TITLE) > 67 else 16\n\nchart = (\n    alt.layer(\n        contour_layer,\n        grid_circles,\n        cross_lines,\n        great_circles,\n        prim_circle,\n        tick_marks,\n        dir_labels,\n        poles_layer,\n        mean_markers,\n        mean_labels,\n    )\n    .properties(\n        width=540,\n        height=540,\n        background=PAGE_BG,\n        padding={\"left\": 6, \"right\": 6, \"top\": 6, \"bottom\": 6},\n        title=alt.Title(\n            text=TITLE,\n            subtitle=\"Lower-Hemisphere Equal-Area (Schmidt) Projection\",\n            fontSize=title_fontsize,\n            subtitleFontSize=11,\n            color=INK,\n            subtitleColor=INK_SOFT,\n        ),\n    )\n    .configure_view(fill=PAGE_BG, stroke=None)\n)\n\n# Save — square target 2400 × 2400. vl-convert pads outside width/height, so\n# pad (never crop) the saved PNG up to the exact canonical canvas.\nchart.save(f\"plot-{THEME}.png\", scale_factor=4.0)\nchart.save(f\"plot-{THEME}.html\")\n\nfrom PIL import Image\n\n\nTW, TH = 2400, 2400\n_img = Image.open(f\"plot-{THEME}.png\").convert(\"RGB\")\n_w, _h = _img.size\nif _w > TW or _h > TH:\n    raise SystemExit(\n        f\"altair vl-convert produced {_w}×{_h}, exceeds target {TW}×{TH}. \"\n        f\"Shrink chart .properties(width=, height=) values and re-render.\"\n    )\nif _w < TW or _h < TH:\n    _canvas = Image.new(\"RGB\", (TW, TH), PAGE_BG)\n    _canvas.paste(_img, ((TW - _w) // 2, (TH - _h) // 2))\n    _canvas.save(f\"plot-{THEME}.png\")\n"}