{"spec_id":"ternary-density","library":"plotnine","language":"python","code":"\"\"\" anyplot.ai\nternary-density: Ternary Density Plot\nLibrary: plotnine 0.15.4 | Python 3.13.13\nQuality: 80/100 | Updated: 2026-05-19\n\"\"\"\n\nimport numpy as np\nimport pandas as pd\nfrom plotnine import (\n    aes,\n    coord_fixed,\n    element_blank,\n    element_text,\n    geom_density_2d,\n    geom_polygon,\n    geom_segment,\n    geom_text,\n    geom_tile,\n    ggplot,\n    labs,\n    scale_fill_cmap,\n    theme,\n    theme_void,\n)\nfrom scipy.stats import gaussian_kde\n\n\n# Generate synthetic compositional data (sediment: sand/silt/clay)\nnp.random.seed(42)\n\n# Create three clusters of compositions\n# Cluster 1: Sandy sediments (high sand)\nn1 = 300\nsand1 = np.random.beta(5, 2, n1) * 60 + 35\nsilt1 = np.random.beta(2, 3, n1) * (100 - sand1) * 0.6\nclay1 = 100 - sand1 - silt1\n\n# Cluster 2: Silty sediments (high silt)\nn2 = 250\nsilt2 = np.random.beta(5, 2, n2) * 50 + 40\nsand2 = np.random.beta(2, 4, n2) * (100 - silt2) * 0.5\nclay2 = 100 - sand2 - silt2\n\n# Cluster 3: Clay-rich sediments (high clay)\nn3 = 250\nclay3 = np.random.beta(4, 2, n3) * 45 + 40\nsand3 = np.random.beta(2, 5, n3) * (100 - clay3) * 0.4\nsilt3 = 100 - sand3 - clay3\n\n# Combine all clusters\nsand = np.concatenate([sand1, sand2, sand3])\nsilt = np.concatenate([silt1, silt2, silt3])\nclay = np.concatenate([clay1, clay2, clay3])\n\n# Ensure non-negative values and normalize to sum to 100\nsand = np.clip(sand, 0, 100)\nsilt = np.clip(silt, 0, 100)\nclay = np.clip(clay, 0, 100)\ntotal = sand + silt + clay\nsand = sand / total * 100\nsilt = silt / total * 100\nclay = clay / total * 100\n\n# Convert ternary to Cartesian coordinates\n# Formula: x = 0.5 * (2*b + c) / total, y = (sqrt(3)/2) * c / total\n# Where a=sand (bottom-left), b=silt (bottom-right), c=clay (top)\nx_data = 0.5 * (2 * silt + clay) / 100\ny_data = (np.sqrt(3) / 2) * clay / 100\n\n# Compute KDE on the data points\npoints = np.vstack([x_data, y_data])\nkde = gaussian_kde(points, bw_method=\"silverman\")\n\n# Create grid for density estimation\ngrid_resolution = 150\nx_grid = np.linspace(0, 1, grid_resolution)\ny_grid = np.linspace(0, np.sqrt(3) / 2, grid_resolution)\nxx, yy = np.meshgrid(x_grid, y_grid)\n\n# Evaluate KDE on grid\npositions = np.vstack([xx.ravel(), yy.ravel()])\ndensity = kde(positions).reshape(xx.shape)\n\n# Create mask for valid ternary region (inside triangle)\n# Convert grid points back to ternary to check bounds\n# c = y * 2/sqrt(3), b = x - c/2, a = 1 - b - c\ncc = yy * 2 / np.sqrt(3)\nbb = xx - cc / 2\naa = 1 - bb - cc\nmask = (aa >= -0.001) & (bb >= -0.001) & (cc >= -0.001) & (aa <= 1.001) & (bb <= 1.001) & (cc <= 1.001)\n\n# Apply mask (set values outside triangle to NaN)\ndensity_masked = np.where(mask, density, np.nan)\n\n# Create DataFrame for density grid\ndensity_df = pd.DataFrame({\"x\": xx.ravel(), \"y\": yy.ravel(), \"density\": density_masked.ravel()}).dropna()\n\n# Triangle vertices for frame\nvertices = pd.DataFrame({\"x\": [0, 1, 0.5, 0], \"y\": [0, 0, np.sqrt(3) / 2, 0]})\n\n# Grid lines at 20% intervals\ngrid_lines = []\nfor pct in [0.2, 0.4, 0.6, 0.8]:\n    # Lines parallel to bottom (constant clay)\n    x1 = 0.5 * (2 * 0 + pct)\n    y1 = (np.sqrt(3) / 2) * pct\n    x2 = 0.5 * (2 * (1 - pct) + pct)\n    y2 = (np.sqrt(3) / 2) * pct\n    grid_lines.append({\"x\": x1, \"y\": y1, \"xend\": x2, \"yend\": y2})\n\n    # Lines parallel to left side (constant silt)\n    x1 = 0.5 * (2 * pct + (1 - pct))\n    y1 = (np.sqrt(3) / 2) * (1 - pct)\n    x2 = 0.5 * (2 * pct + 0)\n    y2 = 0\n    grid_lines.append({\"x\": x1, \"y\": y1, \"xend\": x2, \"yend\": y2})\n\n    # Lines parallel to right side (constant sand)\n    x1 = 0.5 * (2 * 0 + (1 - pct))\n    y1 = (np.sqrt(3) / 2) * (1 - pct)\n    x2 = 0.5 * (2 * (1 - pct) + 0)\n    y2 = 0\n    grid_lines.append({\"x\": x1, \"y\": y1, \"xend\": x2, \"yend\": y2})\n\ngrid_df = pd.DataFrame(grid_lines)\n\n# Tick labels along edges\ntick_labels = []\nlabel_offset = 0.06\nfor pct in [0, 20, 40, 60, 80, 100]:\n    frac = pct / 100\n    # Sand axis (left edge going up)\n    x = 0.5 * (2 * 0 + frac)\n    y = (np.sqrt(3) / 2) * frac\n    tick_labels.append({\"x\": x - label_offset, \"y\": y, \"label\": str(pct)})\n\n    # Silt axis (bottom edge)\n    x = 0.5 * (2 * frac + 0)\n    y = 0\n    tick_labels.append({\"x\": x, \"y\": y - label_offset * 0.8, \"label\": str(pct)})\n\n    # Clay axis (right edge going up)\n    x = 0.5 * (2 * (1 - frac) + frac)\n    y = (np.sqrt(3) / 2) * frac\n    tick_labels.append({\"x\": x + label_offset, \"y\": y, \"label\": str(pct)})\n\ntick_df = pd.DataFrame(tick_labels)\n\n# Vertex labels\nvertex_labels = pd.DataFrame(\n    {\n        \"x\": [0 - 0.02, 1 + 0.02, 0.5],\n        \"y\": [0 - 0.1, 0 - 0.1, np.sqrt(3) / 2 + 0.08],\n        \"label\": [\"Sand (%)\", \"Silt (%)\", \"Clay (%)\"],\n    }\n)\n\n# Build the plot\n# Calculate tile size for proper coverage\ntile_size = 1.0 / grid_resolution\n\nplot = (\n    ggplot()\n    # Density heatmap using tiles\n    + geom_tile(\n        data=density_df,\n        mapping=aes(x=\"x\", y=\"y\", fill=\"density\"),\n        width=tile_size * 1.1,\n        height=tile_size * 1.1,\n        alpha=0.85,\n    )\n    + scale_fill_cmap(cmap_name=\"viridis\", name=\"Density\")\n    # Contour lines for key density levels using geom_density_2d\n    + geom_density_2d(\n        data=pd.DataFrame({\"x\": x_data, \"y\": y_data}),\n        mapping=aes(x=\"x\", y=\"y\"),\n        color=\"white\",\n        size=0.8,\n        alpha=0.7,\n        levels=6,\n    )\n    # Grid lines beneath density (already visible due to alpha)\n    + geom_segment(\n        data=grid_df, mapping=aes(x=\"x\", y=\"y\", xend=\"xend\", yend=\"yend\"), color=\"#ffffff\", size=0.8, alpha=0.3\n    )\n    # Triangle frame on top\n    + geom_polygon(data=vertices, mapping=aes(x=\"x\", y=\"y\"), fill=None, color=\"#306998\", size=2.5)\n    # Tick labels\n    + geom_text(data=tick_df, mapping=aes(x=\"x\", y=\"y\", label=\"label\"), size=14, color=\"#666666\")\n    # Vertex labels\n    + geom_text(\n        data=vertex_labels, mapping=aes(x=\"x\", y=\"y\", label=\"label\"), size=18, fontweight=\"bold\", color=\"#306998\"\n    )\n    # Title and theme\n    + labs(title=\"ternary-density · plotnine · pyplots.ai\")\n    + coord_fixed(ratio=1)\n    + theme_void()\n    + theme(\n        figure_size=(12, 12),\n        plot_title=element_text(size=28, ha=\"center\", weight=\"bold\"),\n        plot_margin=0.02,\n        legend_position=\"right\",\n        legend_title=element_text(size=18),\n        legend_text=element_text(size=14),\n        legend_key_size=40,\n        legend_background=element_blank(),\n    )\n)\n\n# Save\nplot.save(\"plot.png\", dpi=300, verbose=False)\n"}