{"spec_id":"ternary-density","library":"letsplot","language":"python","code":"\"\"\" anyplot.ai\nternary-density: Ternary Density Plot\nLibrary: letsplot 4.9.0 | Python 3.13.13\nQuality: 87/100 | Updated: 2026-05-19\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_path,\n    geom_point,\n    geom_polygon,\n    geom_segment,\n    geom_text,\n    ggplot,\n    ggsave,\n    ggsize,\n    labs,\n    layer_tooltips,\n    scale_fill_viridis,\n    theme,\n)\n\n\nLetsPlot.setup_html()\n\n# Theme 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\"\n\n# Data — synthetic compositional data (sediment: sand/silt/clay)\nnp.random.seed(42)\n\n# Three clusters via Dirichlet distribution\nalpha1 = np.array([8, 2, 1])\ncomp1 = np.random.dirichlet(alpha1, 180) * 100\n\nalpha2 = np.array([2, 7, 2])\ncomp2 = np.random.dirichlet(alpha2, 160) * 100\n\nalpha3 = np.array([1, 2, 8])\ncomp3 = np.random.dirichlet(alpha3, 160) * 100\n\ncompositions = np.vstack([comp1, comp2, comp3])\nsand = compositions[:, 0]\nsilt = compositions[:, 1]\nclay = compositions[:, 2]\n\n# Convert ternary to Cartesian coordinates\n# bottom-left = Sand, bottom-right = Silt, top = Clay\ntotal = sand + silt + clay\nb_norm = silt / total\nc_norm = clay / total\n\nx_data = 0.5 * (2 * b_norm + c_norm)\ny_data = (np.sqrt(3) / 2) * c_norm\n\n# Density grid\ngrid_res = 100\nx_grid = np.linspace(0, 1, grid_res)\ny_grid = np.linspace(0, np.sqrt(3) / 2, grid_res)\nX, Y = np.meshgrid(x_grid, y_grid)\n\n# 2D Gaussian KDE (Scott's rule)\nn = len(x_data)\nbw = n ** (-1.0 / 6)\nbw_x = np.std(x_data) * bw\nbw_y = np.std(y_data) * bw\n\nZ = np.zeros_like(X)\nfor i in range(n):\n    dx = (X - x_data[i]) / bw_x\n    dy = (Y - y_data[i]) / bw_y\n    Z += np.exp(-0.5 * (dx**2 + dy**2))\nZ /= n * 2 * np.pi * bw_x * bw_y\n\n# Mask points outside the equilateral triangle\nsqrt3 = np.sqrt(3)\nmask = (Y >= 0) & (Y <= sqrt3 * X + 1e-6) & (Y <= sqrt3 * (1 - X) + 1e-6)\n\n# Density polygons dataframe\npolygon_data = []\npoly_id = 0\ndx = x_grid[1] - x_grid[0]\ndy = y_grid[1] - y_grid[0]\noverlap = 1.05\n\nfor i in range(grid_res):\n    for j in range(grid_res):\n        if mask[i, j] and Z[i, j] > 0:\n            cx, cy = X[i, j], Y[i, j]\n            hdx = dx * overlap / 2\n            hdy = dy * overlap / 2\n            corners_x = [cx - hdx, cx + hdx, cx + hdx, cx - hdx, cx - hdx]\n            corners_y = [cy - hdy, cy - hdy, cy + hdy, cy + hdy, cy - hdy]\n            for k in range(5):\n                polygon_data.append({\"x\": corners_x[k], \"y\": corners_y[k], \"density\": Z[i, j], \"id\": poly_id})\n            poly_id += 1\n\ndf_polygons = pd.DataFrame(polygon_data)\n\n# Contour lines via marching squares at 25%, 50%, 75% density levels\nz_masked = Z.copy()\nz_masked[~mask] = 0\nz_min, z_max = z_masked[mask].min(), z_masked[mask].max()\ncontour_levels = [z_min + (z_max - z_min) * p for p in [0.25, 0.5, 0.75]]\n\ncontour_data = []\nfor level in contour_levels:\n    for i in range(grid_res - 1):\n        for j in range(grid_res - 1):\n            corners = [Z[i, j], Z[i, j + 1], Z[i + 1, j + 1], Z[i + 1, j]]\n            corners_mask = [mask[i, j], mask[i, j + 1], mask[i + 1, j + 1], mask[i + 1, j]]\n            if not all(corners_mask):\n                continue\n            above = [c >= level for c in corners]\n            if all(above) or not any(above):\n                continue\n            x0, x1 = x_grid[j], x_grid[j + 1]\n            y0, y1 = y_grid[i], y_grid[i + 1]\n            pts = []\n            if above[0] != above[1]:\n                t = (level - corners[0]) / (corners[1] - corners[0] + 1e-10)\n                pts.append((x0 + t * (x1 - x0), y0))\n            if above[1] != above[2]:\n                t = (level - corners[1]) / (corners[2] - corners[1] + 1e-10)\n                pts.append((x1, y0 + t * (y1 - y0)))\n            if above[2] != above[3]:\n                t = (level - corners[2]) / (corners[3] - corners[2] + 1e-10)\n                pts.append((x1 - t * (x1 - x0), y1))\n            if above[3] != above[0]:\n                t = (level - corners[3]) / (corners[0] - corners[3] + 1e-10)\n                pts.append((x0, y1 - t * (y1 - y0)))\n            if len(pts) == 2:\n                contour_data.append({\"x\": pts[0][0], \"y\": pts[0][1], \"xend\": pts[1][0], \"yend\": pts[1][1]})\n\ndf_contours = pd.DataFrame(contour_data) if contour_data else pd.DataFrame(columns=[\"x\", \"y\", \"xend\", \"yend\"])\n\n# Triangle outline\ntri_x = [0, 1, 0.5, 0]\ntri_y = [0, 0, sqrt3 / 2, 0]\ndf_triangle = pd.DataFrame({\"x\": tri_x, \"y\": tri_y})\n\n# Grid lines inside the triangle\ngrid_lines = []\nfor pct in [0.2, 0.4, 0.6, 0.8]:\n    # Parallel to bottom edge (constant clay %)\n    y_line = pct * sqrt3 / 2\n    grid_lines.append({\"x\": pct / 2, \"xend\": 1 - pct / 2, \"y\": y_line, \"yend\": y_line})\n    # Parallel to left edge (constant silt %)\n    grid_lines.append({\"x\": pct, \"xend\": 1 - 0.5 * pct, \"y\": 0.0, \"yend\": pct * sqrt3 / 2})\n    # Parallel to right edge (constant sand %)\n    grid_lines.append({\"x\": 1 - pct, \"xend\": 0.5 * pct, \"y\": 0.0, \"yend\": pct * sqrt3 / 2})\n\ndf_grid = pd.DataFrame(grid_lines)\n\n# Vertex labels\nlabels_data = pd.DataFrame(\n    {\"x\": [-0.06, 1.06, 0.5], \"y\": [-0.05, -0.05, sqrt3 / 2 + 0.06], \"label\": [\"Sand\", \"Silt\", \"Clay\"]}\n)\n\n# Scatter data for interactive HTML tooltips showing composition at each point\ndf_scatter = pd.DataFrame(\n    {\"x\": x_data, \"y\": y_data, \"Sand\": sand.round(1), \"Silt\": silt.round(1), \"Clay\": clay.round(1)}\n)\n\n# Plot\nplot = (\n    ggplot()\n    + geom_polygon(aes(x=\"x\", y=\"y\", fill=\"density\", group=\"id\"), data=df_polygons, color=None, alpha=0.9)\n    + scale_fill_viridis(name=\"KDE Density\", option=\"viridis\")\n    + geom_segment(aes(x=\"x\", y=\"y\", xend=\"xend\", yend=\"yend\"), data=df_grid, color=INK_SOFT, size=1.0, alpha=0.5)\n    + (\n        geom_segment(\n            aes(x=\"x\", y=\"y\", xend=\"xend\", yend=\"yend\"),\n            data=df_contours,\n            color=\"white\",\n            size=2.5,\n            alpha=0.9,\n            linetype=\"dashed\",\n        )\n        if len(df_contours) > 0\n        else geom_path(aes(x=\"x\", y=\"y\"), data=pd.DataFrame({\"x\": [], \"y\": []}))\n    )\n    + geom_path(aes(x=\"x\", y=\"y\"), data=df_triangle, color=INK, size=2.0)\n    + geom_point(\n        aes(x=\"x\", y=\"y\"),\n        data=df_scatter,\n        color=\"white\",\n        size=1,\n        alpha=0.01,\n        tooltips=layer_tooltips().line(\"Sand: @Sand%\").line(\"Silt: @Silt%\").line(\"Clay: @Clay%\"),\n    )\n    + geom_text(aes(x=\"x\", y=\"y\", label=\"label\"), data=labels_data, color=INK, size=14, fontface=\"bold\")\n    + labs(title=\"Sediment Composition · ternary-density · python · letsplot · anyplot.ai\", x=\"\", y=\"\")\n    + coord_fixed(ratio=1)\n    + theme(\n        plot_background=element_rect(fill=PAGE_BG, color=PAGE_BG),\n        panel_background=element_rect(fill=PAGE_BG),\n        panel_grid=element_blank(),\n        axis_title=element_blank(),\n        axis_text=element_blank(),\n        axis_ticks=element_blank(),\n        axis_line=element_blank(),\n        plot_title=element_text(size=24, face=\"bold\", color=INK),\n        legend_background=element_rect(fill=ELEVATED_BG, color=INK_SOFT),\n        legend_title=element_text(size=16, color=INK),\n        legend_text=element_text(size=14, color=INK_SOFT),\n    )\n    + ggsize(1600, 900)\n)\n\n# Save\nggsave(plot, f\"plot-{THEME}.png\", path=\".\", scale=3)\nggsave(plot, f\"plot-{THEME}.html\", path=\".\")\n"}