{"spec_id":"stereonet-equal-area","library":"plotnine","language":"python","code":"\"\"\" anyplot.ai\nstereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\nLibrary: plotnine 0.15.7 | Python 3.13.13\nQuality: 90/100 | Updated: 2026-06-16\n\"\"\"\n\nimport os\n\nimport numpy as np\nimport pandas as pd\nfrom plotnine import (\n    aes,\n    after_stat,\n    coord_fixed,\n    element_blank,\n    element_rect,\n    element_text,\n    geom_density_2d,\n    geom_path,\n    geom_point,\n    geom_segment,\n    geom_text,\n    ggplot,\n    labs,\n    scale_alpha_continuous,\n    scale_color_manual,\n    scale_shape_manual,\n    scale_x_continuous,\n    scale_y_continuous,\n    theme,\n)\n\n\n# Theme tokens (see prompts/default-style-guide.md \"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\"\nINK_MUTED = \"#6B6A63\" if THEME == \"light\" else \"#A8A79F\"\n\n# Imprint palette — Bedding=brand green (always first), Fault=lavender, Joint=blue\nIMPRINT = {\"Bedding\": \"#009E73\", \"Fault\": \"#C475FD\", \"Joint\": \"#4467A3\"}\nSHAPES = {\"Bedding\": \"o\", \"Fault\": \"D\", \"Joint\": \"s\"}\n\n# Data — geological field measurements (strike/dip for bedding, faults, joints)\nnp.random.seed(42)\n\nbedding_strike = np.random.normal(45, 12, 40)\nbedding_dip = np.clip(np.random.normal(30, 8, 40), 5, 85)\n\nfault_strike = np.random.normal(150, 15, 25)\nfault_dip = np.clip(np.random.normal(65, 10, 25), 10, 88)\n\njoint_strike = np.random.normal(280, 20, 35)\njoint_dip = np.clip(np.random.normal(75, 8, 35), 15, 89)\n\nstrikes = np.concatenate([bedding_strike, fault_strike, joint_strike]) % 360\ndips = np.concatenate([bedding_dip, fault_dip, joint_dip])\nfeature_types = [\"Bedding\"] * 40 + [\"Fault\"] * 25 + [\"Joint\"] * 35\n\n# Primitive circle radius for the Schmidt equal-area projection\nr_prim = np.sqrt(2)\n\n# Compute poles to planes (equal-area lower-hemisphere projection)\npole_trend = np.radians((strikes + 90) % 360)\npole_plunge = np.radians(90 - dips)\npole_r = np.sqrt(2) * np.sin((np.pi / 2 - pole_plunge) / 2)\npole_x = pole_r * np.sin(pole_trend)\npole_y = pole_r * np.cos(pole_trend)\npoles_df = pd.DataFrame({\"x\": pole_x, \"y\": pole_y, \"feature_type\": feature_types})\n\n# Compute great circles for a few representative planes per feature type\ngc_rows = []\ngc_indices = {\"Bedding\": [0, 10, 20, 30], \"Fault\": [40, 48, 56], \"Joint\": [65, 75, 85]}\ngc_id = 0\nfor ftype, indices in gc_indices.items():\n    for idx in indices:\n        if idx >= len(strikes):\n            continue\n        strike_rad = np.radians(strikes[idx])\n        dip_rad = np.radians(dips[idx])\n        strike_vec = np.array([np.sin(strike_rad), np.cos(strike_rad), 0.0])\n        dip_dir_rad = strike_rad + np.pi / 2\n        dip_vec = np.array(\n            [np.sin(dip_dir_rad) * np.cos(dip_rad), np.cos(dip_dir_rad) * np.cos(dip_rad), -np.sin(dip_rad)]\n        )\n        for a in np.linspace(-np.pi / 2, np.pi / 2, 181):\n            pt = np.cos(a) * dip_vec + np.sin(a) * strike_vec\n            if pt[2] > 0:\n                pt = -pt\n            horiz = np.sqrt(pt[0] ** 2 + pt[1] ** 2)\n            plunge = np.arctan2(-pt[2], horiz)\n            trend = np.arctan2(pt[0], pt[1])\n            r = np.sqrt(2) * np.sin((np.pi / 2 - plunge) / 2)\n            x, y = r * np.sin(trend), r * np.cos(trend)\n            if x**2 + y**2 <= r_prim**2 * 1.01:\n                gc_rows.append({\"x\": x, \"y\": y, \"feature_type\": ftype, \"gc_id\": f\"{ftype}_{gc_id}\"})\n        gc_id += 1\n\ngc_df = pd.DataFrame(gc_rows)\n\n# Stereonet net — primitive circle\ncircle_angles = np.linspace(0, 2 * np.pi, 361)\nprim_df = pd.DataFrame({\"x\": r_prim * np.cos(circle_angles), \"y\": r_prim * np.sin(circle_angles)})\n\n# Equal-area net grid — small circles at 30 deg dip intervals\ngrid_rows = []\nfor dip_interval in range(30, 90, 30):\n    plunge_rad = np.radians(90 - dip_interval)\n    r_circle = np.sqrt(2) * np.sin((np.pi / 2 - plunge_rad) / 2)\n    for angle in np.linspace(0, 2 * np.pi, 181):\n        grid_rows.append(\n            {\"x\": r_circle * np.cos(angle), \"y\": r_circle * np.sin(angle), \"grid_id\": f\"dip_{dip_interval}\"}\n        )\n\n# Radial lines at 30 deg azimuth intervals\nfor az in range(0, 360, 30):\n    az_rad = np.radians(az)\n    for t in np.linspace(0, r_prim, 50):\n        grid_rows.append({\"x\": t * np.sin(az_rad), \"y\": t * np.cos(az_rad), \"grid_id\": f\"az_{az}\"})\n\ngrid_df = pd.DataFrame(grid_rows)\n\n# Degree tick marks every 10 degrees around the perimeter\ntick_rows = []\nfor deg in range(0, 360, 10):\n    rad = np.radians(deg)\n    inner, outer = r_prim * 0.97, r_prim * 1.0\n    tick_rows.append(\n        {\"x1\": inner * np.sin(rad), \"y1\": inner * np.cos(rad), \"x2\": outer * np.sin(rad), \"y2\": outer * np.cos(rad)}\n    )\ntick_df = pd.DataFrame(tick_rows)\n\n# Cardinal direction labels\ndir_labels = []\nfor deg, label in [(0, \"N\"), (90, \"E\"), (180, \"S\"), (270, \"W\")]:\n    rad = np.radians(deg)\n    offset = r_prim * 1.13\n    dir_labels.append({\"x\": offset * np.sin(rad), \"y\": offset * np.cos(rad), \"label\": label})\ndir_df = pd.DataFrame(dir_labels)\n\n# Degree labels every 30 degrees (excluding cardinal directions)\ndeg_labels = []\nfor deg in range(0, 360, 30):\n    if deg in (0, 90, 180, 270):\n        continue\n    rad = np.radians(deg)\n    offset = r_prim * 1.09\n    deg_labels.append({\"x\": offset * np.sin(rad), \"y\": offset * np.cos(rad), \"label\": f\"{deg}°\"})\ndeg_label_df = pd.DataFrame(deg_labels)\n\n# Mean strike/dip annotation per cluster (geological context)\nannotations = []\nfor ftype in [\"Bedding\", \"Fault\", \"Joint\"]:\n    mask = poles_df[\"feature_type\"] == ftype\n    cx, cy = poles_df.loc[mask, \"x\"].mean(), poles_df.loc[mask, \"y\"].mean()\n    mean_strike = strikes[np.array(feature_types) == ftype].mean()\n    mean_dip = dips[np.array(feature_types) == ftype].mean()\n    annotations.append({\"x\": cx, \"y\": cy - 0.14, \"feature_type\": ftype, \"label\": f\"{mean_strike:.0f}°/{mean_dip:.0f}°\"})\nannot_df = pd.DataFrame(annotations)\n\n# Plot — layered grammar of graphics on the equal-area net\nplot = (\n    ggplot()\n    # Equal-area net grid (subtle, theme-adaptive)\n    + geom_path(aes(x=\"x\", y=\"y\", group=\"grid_id\"), data=grid_df, color=INK, size=0.3, alpha=0.15)\n    # Primitive circle (projection boundary)\n    + geom_path(aes(x=\"x\", y=\"y\"), data=prim_df, color=INK, size=1.2)\n    # Perimeter degree ticks\n    + geom_segment(aes(x=\"x1\", y=\"y1\", xend=\"x2\", yend=\"y2\"), data=tick_df, color=INK_SOFT, size=0.6)\n    # Pole concentration density contours\n    + geom_density_2d(\n        aes(x=\"x\", y=\"y\", alpha=after_stat(\"level\")),\n        data=poles_df,\n        color=INK_MUTED,\n        size=0.6,\n        linetype=\"dashed\",\n        show_legend=False,\n    )\n    + scale_alpha_continuous(range=(0.25, 0.7))\n    # Great circles for representative planes\n    + geom_path(aes(x=\"x\", y=\"y\", color=\"feature_type\", group=\"gc_id\"), data=gc_df, size=0.9, alpha=0.65)\n    # Poles to planes\n    + geom_point(\n        aes(x=\"x\", y=\"y\", color=\"feature_type\", shape=\"feature_type\"), data=poles_df, size=3.5, alpha=0.9, stroke=0.5\n    )\n    # Cardinal direction labels\n    + geom_text(aes(x=\"x\", y=\"y\", label=\"label\"), data=dir_df, size=7, fontweight=\"bold\", color=INK)\n    # Perimeter degree labels\n    + geom_text(aes(x=\"x\", y=\"y\", label=\"label\"), data=deg_label_df, size=3.6, color=INK_SOFT)\n    # Mean strike/dip per cluster\n    + geom_text(\n        aes(x=\"x\", y=\"y\", label=\"label\", color=\"feature_type\"),\n        data=annot_df,\n        size=4.0,\n        fontstyle=\"italic\",\n        fontweight=\"bold\",\n        show_legend=False,\n    )\n    + scale_color_manual(name=\"Feature type\", values=IMPRINT)\n    + scale_shape_manual(name=\"Feature type\", values=SHAPES)\n    + coord_fixed(ratio=1)\n    + scale_x_continuous(limits=(-1.85, 1.85))\n    + scale_y_continuous(limits=(-1.85, 1.85))\n    + labs(title=\"stereonet-equal-area · python · plotnine · anyplot.ai\")\n    + theme(\n        figure_size=(6, 6),\n        plot_background=element_rect(fill=PAGE_BG, color=PAGE_BG),\n        panel_background=element_rect(fill=PAGE_BG, color=PAGE_BG),\n        legend_background=element_rect(fill=ELEVATED_BG, color=INK_SOFT),\n        legend_key=element_rect(fill=ELEVATED_BG, color=ELEVATED_BG),\n        plot_title=element_text(size=13, ha=\"center\", color=INK),\n        legend_title=element_text(size=11, weight=\"bold\", color=INK),\n        legend_text=element_text(size=10, color=INK_SOFT),\n        legend_position=\"bottom\",\n        axis_title=element_blank(),\n        axis_text=element_blank(),\n        axis_ticks=element_blank(),\n        axis_line=element_blank(),\n        panel_grid_major=element_blank(),\n        panel_grid_minor=element_blank(),\n        panel_border=element_blank(),\n    )\n)\n\n# Save (2400 x 2400 px square at dpi=400)\nplot.save(f\"plot-{THEME}.png\", dpi=400, width=6, height=6, units=\"in\")\n"}