Skip to content

Instantly share code, notes, and snippets.

@rutj3

rutj3/soil_texture_warped.ipynb

Last active April 25, 2025 13:22
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save rutj3/1dcae00e2ac56f392b29dc12fd74f480 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save rutj3/1dcae00e2ac56f392b29dc12fd74f480 to your computer and use it in GitHub Desktop.
Download ZIP
Soil texture warped
Display the source blob
Display the rendered blob
Raw
soil_texture_warped.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "806ea258-de01-4bd0-a733-39178ea7c88c",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import log, exp, sqrt, cos, sin\n",
"import numba\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"from collections import defaultdict\n",
"import matplotlib.patheffects as PathEffects\n",
"plt.style.use(\"seaborn-v0_8\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "58b46be3-c307-44f1-b35c-fd18890581f9",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def get_geometric_props(fy, ft, fd):\n",
"\n",
" # # Shirazi & Boersma (1984) in mm\n",
" # clay_lo = 0\n",
" # clay_hi = 0.002\n",
" # silt_lo = 0.002\n",
" # silt_hi = 0.05\n",
" # sand_lo = 0.05\n",
" # sand_hi = 2.00\n",
" \n",
" # clay_avg = (clay_lo + clay_hi) / 2\n",
" # silt_avg = (silt_lo + silt_hi) / 2\n",
" # sand_avg = (sand_lo + sand_hi) / 2\n",
"\n",
" # Boersma & Hart (1988) in μm\n",
" clay_lo = 0.01\n",
" clay_hi = 2.0\n",
" silt_lo = 2.0\n",
" silt_hi = 50.0\n",
" sand_lo = 50.0\n",
" sand_hi = 2000.0\n",
" \n",
" clay_avg = sqrt(clay_lo * clay_hi)\n",
" silt_avg = sqrt(silt_lo * silt_hi)\n",
" sand_avg = sqrt(sand_lo * sand_hi)\n",
" \n",
" a = 0.01 * (fy * log(clay_avg) + ft * log(silt_avg) + fd * log(sand_avg))\n",
" b = sqrt(np.maximum((0.01 * (fy * log(clay_avg)**2 + ft * log(silt_avg)**2 + fd * log(sand_avg)**2)) - a**2, 0))\n",
" dg = exp(a)\n",
" σg = exp(b)\n",
"\n",
" return dg, σg"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3446bc73-b1aa-4332-a3e2-522aa9b4e294",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"clay_boundaries = []\n",
"silt_boundaries = []\n",
"sand_boundaries = []\n",
"\n",
"line_step = 10\n",
"line_res = .01\n",
"\n",
"for i, fy in enumerate(range(0,100+line_step,line_step)):\n",
" \n",
" ft = np.arange(0, 100+line_res-fy, line_res)\n",
" fy = np.full_like(ft, fy)\n",
" fd = 100 - fy - ft\n",
" dg, σg = get_geometric_props(fy, ft, fd)\n",
" clay_boundaries.append((dg, σg))\n",
" \n",
"for i, ft in enumerate(range(0,100+line_step,line_step)):\n",
" \n",
" fy = np.arange(0, 100+line_res-ft, line_res)\n",
" ft = np.full_like(fy, ft)\n",
" fd = 100 - fy - ft\n",
" dg, σg = get_geometric_props(fy, ft, fd)\n",
" silt_boundaries.append((dg, σg))\n",
" \n",
"for i, fd in enumerate(range(0,100+line_step,line_step)):\n",
" \n",
" fy = np.arange(0,100+line_res-fd, line_res)\n",
" fd = np.full_like(fy, fd)\n",
" ft = 100 - fy - fd\n",
" \n",
" dg, σg = get_geometric_props(fy, ft, fd)\n",
" sand_boundaries.append((dg, σg))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "831a47c6-6cab-47b6-beac-0984fe4f25ad",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# name, (fy, ft, fd)\n",
"class_nodes = { \n",
" \"b\": (0, 0, 100),\n",
" \"d\": (0, 30, 70),\n",
" \"f\": (0, 50, 50),\n",
" \"m\": (0, 80, 20),\n",
" \"q\": (0, 100, 0),\n",
" \"c\": (10, 0, 90),\n",
" \"g\": (20, 0, 80),\n",
" \"v\": (40, 40, 20),\n",
" \"w\": (40, 60, 0),\n",
" \"y\": (60, 40, 0),\n",
" \"z\": (100, 0, 0),\n",
" \"a\": (0, 15, 85),\n",
" \"e\": (15, 0, 85),\n",
" \"h\": (20, 28, 52),\n",
" \"i\": (7, 41, 52),\n",
" \"j\": (7, 50, 43),\n",
" \"k\": (27, 28, 45),\n",
" \"l\": (27, 50, 23),\n",
" \"n\": (27, 73, 0),\n",
" \"o\": (12, 88, 0),\n",
" \"p\": (12, 80, 8),\n",
" \"r\": (35, 0, 65),\n",
" \"s\": (35, 20, 45),\n",
" \"t\": (27, 53, 20),\n",
" \"u\": (40, 15, 45),\n",
" \"x\": (55, 0, 45), \n",
"}\n",
"\n",
"class_polys = {\n",
" \"1) sand\": [\"b\",\"c\",\"a\"],\n",
" \"2) loamy sand\": [\"a\",\"c\",\"e\",\"d\"],\n",
" \"3) sandy loam\": [\"d\",\"e\",\"g\",\"h\",\"i\",\"j\",\"f\"],\n",
" \"4) loam\": [\"j\",\"i\",\"h\",\"k\",\"l\"],\n",
" \"5) silt loam\": [\"f\",\"j\",\"l\",\"n\",\"o\",\"p\",\"m\"],\n",
" \"6) silt\": [\"m\",\"p\",\"o\",\"q\"],\n",
" \"7) sandy clay loam\": [\"g\",\"r\",\"s\",\"k\",\"h\"],\n",
" \"8) clay loam\": [\"k\",\"u\",\"v\",\"t\",\"l\"],\n",
" \"9) silt clay loam\": [\"t\",\"v\",\"w\",\"n\"],\n",
" \"10) sandy clay\": [\"r\",\"x\",\"u\",\"s\"],\n",
" \"11) silty clay\": [\"v\",\"y\",\"w\"],\n",
" \"12) clay\": [\"x\",\"z\",\"y\",\"v\",\"u\"],\n",
"}\n",
"\n",
"# clay/silt/sand\n",
"class_label_pos = {\n",
" \"1) sand\": (4, 4),\n",
" \"2) loamy sand\": (5, 13),\n",
" \"3) sandy loam\": (11, 23),\n",
" \"4) loam\": (20, 40),\n",
" \"5) silt loam\": (16, 65),\n",
" \"6) silt\": (5, 87),\n",
" \"7) sandy clay loam\": (26, 15),\n",
" \"8) clay loam\": (33, 33),\n",
" \"9) silt clay loam\": (33, 55),\n",
" \"10) sandy clay\": (41, 8),\n",
" \"11) silty clay\": (47, 47),\n",
" \"12) clay\": (60, 20),\n",
"}\n",
"\n",
"class_polys_full = {}\n",
"\n",
"for name, nodes in class_polys.items():\n",
" class_polys_full[name] = list(map(class_nodes.get, nodes))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ddd2ed5-8b4a-4973-9ce2-8bd9e570d0fe",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"class_bounds = defaultdict(list)\n",
"\n",
"for name, nodes in class_polys_full.items():\n",
" \n",
" for i in range(len(nodes)):\n",
" fy1, ft1, fd1 = nodes[i-1]\n",
" fy2, ft2, fd2 = nodes[i]\n",
" \n",
" distance = sqrt((fy1-fy2)**2 + (ft1-ft2)**2 + (fd1-fd2)**2)\n",
" line_length = int(np.ceil(distance*100))\n",
" \n",
" fy = np.linspace(fy1, fy2, line_length)\n",
" ft = np.linspace(ft1, ft2, line_length)\n",
" fd = np.linspace(fd1, fd2, line_length)\n",
"\n",
" dg, σg = get_geometric_props(fy, ft, fd)\n",
" \n",
" class_bounds[name].append((dg, σg))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc4f92a0-de48-4dbd-a113-019d0df08636",
"metadata": {},
"outputs": [],
"source": [
"axes_poly = [\n",
" ( 0, 0, 99.9),\n",
" (99.9, 0, 0),\n",
" ( 0, 99.9, 0),\n",
" ( 0, 0, 99.9),\n",
"]\n",
"\n",
"n = 1000\n",
"axes_boundary = np.vstack((\n",
" \n",
" np.stack((\n",
" np.linspace(0, 100, n),\n",
" np.linspace(0, 0, n),\n",
" np.linspace(100, 0, n),\n",
" ), axis=1),\n",
" \n",
" np.stack((\n",
" np.linspace(100, 0, n),\n",
" np.linspace(0, 100, n),\n",
" np.linspace(0, 0, n),\n",
" ), axis=1),\n",
" \n",
" np.stack((\n",
" np.linspace( 0, 0, n),\n",
" np.linspace(100, 0, n),\n",
" np.linspace(0, 100, n),\n",
" ), axis=1)\n",
"))\n",
"\n",
"boundary_dg, boundary_σg = get_geometric_props(\n",
" axes_boundary[:, 0],\n",
" axes_boundary[:, 1],\n",
" axes_boundary[:, 2],\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "abfd8222-eff6-4b9e-92f3-c137b8cba970",
"metadata": {},
"outputs": [],
"source": [
"def annotate(\n",
" ax, label_min=20, label_max=80, label_step=20, position_offset=35, position_rotation_offset=0, txt_rotation_offset=0, \n",
" color=\"C0\", plot_marker=True, axis=\"Silt\", txt=\"{value:1.0f}%\", \n",
" coord_conversion=get_geometric_props, marker_angle=None, ha=None, va=None, align_with=\"boundary\", stroke=False,\n",
" debug=False,\n",
"):\n",
"\n",
" data_to_figure = ax.transData.transform\n",
" data_to_axes = lambda x: ax.transAxes.inverted().transform(data_to_figure(x))\n",
"\n",
" for value in range(label_min, label_max+1, label_step):\n",
"\n",
" if align_with == \"boundary\":\n",
" # finite difference along the boundary for angle\n",
" if axis.lower() == \"silt\":\n",
" ft0 = value\n",
" ft1 = max(ft0 - 0.01, 0)\n",
" ft2 = max(ft0 + 0.01, 0)\n",
" \n",
" fd0 = fd1 = fd2 = 0\n",
" \n",
" fy0 = 100 - ft0\n",
" fy1 = 100 - ft1\n",
" fy2 = 100 - ft2\n",
" \n",
" elif axis.lower() == \"clay\":\n",
" fy0 = value\n",
" fy1 = max(fy0 - 0.01, 0)\n",
" fy2 = max(fy0 + 0.01, 0)\n",
" \n",
" ft0 = ft1 = ft2 = 0\n",
" \n",
" fd0 = 100 - fy0\n",
" fd1 = 100 - fy1\n",
" fd2 = 100 - fy2\n",
" \n",
" elif axis.lower() == \"sand\":\n",
" fd0 = value\n",
" fd1 = max(fd0 - 0.01, 0)\n",
" fd2 = max(fd0 + 0.01, 0)\n",
" \n",
" ft0 = 100 - fd0\n",
" ft1 = 100 - fd1\n",
" ft2 = 100 - fd2\n",
" \n",
" fy0 = fy1 = fy2 = 0\n",
"\n",
" elif align_with == \"grid\":\n",
" # finite difference along the gridline for angle\n",
" if axis.lower() == \"silt\":\n",
" ft0 = ft1 = ft2 = value\n",
" \n",
" fd0 = 0\n",
" fd1 = 0.1\n",
" fd2 = 0\n",
" \n",
" fy0 = 100 - ft0\n",
" fy1 = 99.9 - ft1\n",
" fy2 = 100 - ft2\n",
" \n",
" elif axis.lower() == \"clay\":\n",
" fy0 = fy1 = fy2 = value\n",
" \n",
" ft0 = 0\n",
" ft1 = 0\n",
" ft2 = 0.1\n",
" \n",
" fd0 = 100 - fy0\n",
" fd1 = 100 - fy1\n",
" fd2 = 99.9 - fy2\n",
" \n",
" elif axis.lower() == \"sand\":\n",
" fd0 = fd1 = fd2 = value\n",
" \n",
" ft0 = 100. - fd0\n",
" ft1 = 99.9 - fd1\n",
" ft2 = 100. - fd2\n",
" \n",
" fy0 = 0\n",
" fy1 = 0.1\n",
" fy2 = 0\n",
" \n",
" x0, y0 = coord_conversion(fy0, ft0, fd0)\n",
" x1, y1 = coord_conversion(fy1, ft1, fd1)\n",
" x2, y2 = coord_conversion(fy2, ft2, fd2)\n",
" dx = (log(x2) - log(x1) ) / exp(1.4)\n",
" # dx = (x2 - x1)\n",
" dy = (y2 - y1) / 50\n",
"\n",
" ax1, ay1 = data_to_axes((x1, y1))\n",
" ax2, ay2 = data_to_axes((x2, y2))\n",
" # dx = (ax2 - ax1)\n",
" # dy = (ay2 - ay1)\n",
" \n",
" fx1, fy1 = data_to_figure((x1, y1))\n",
" fx2, fy2 = data_to_figure((x2, y2))\n",
" fdx = (fx2 - fx1)\n",
" fdy = (fy2 - fy1)\n",
"\n",
" angle_deg = (np.rad2deg(np.arctan2(-dx, dy))+90) % 360 \n",
" angle_rad = np.deg2rad(angle_deg)\n",
"\n",
" fangle_deg = (np.rad2deg(np.arctan2(-fdx, fdy))+90) % 360\n",
" fangle_rad = np.deg2rad(fangle_deg)\n",
"\n",
" if position_offset == 0:\n",
" sgn = 1\n",
" else:\n",
" sgn = np.sign(position_offset)\n",
" \n",
" label_offset = position_offset\n",
" px0 = (label_offset * cos(angle_rad + np.deg2rad(0) + np.deg2rad(position_rotation_offset)))\n",
" py0 = (label_offset * sin(angle_rad + np.deg2rad(0) + np.deg2rad(position_rotation_offset)))\n",
"\n",
" if plot_marker:\n",
" if marker_angle is None:\n",
" marker_angle = angle_deg + txt_rotation_offset\n",
" \n",
" ax.plot(\n",
" x0, y0, mec=color, ms=15, mfc=\"none\", mew=1., clip_on=False,\n",
" marker=(1, 2, (fangle_deg + txt_rotation_offset + 90*((txt_rotation_offset > 0 and position_offset > 0)*2-1))),\n",
" )\n",
"\n",
" if ha is None:\n",
" ha = \"center\"\n",
" if va is None:\n",
" va = \"center\"\n",
"\n",
" if debug is True:\n",
" bbox_props = dict(facecolor=\"none\", edgecolor=\"m\", lw=1)\n",
" else:\n",
" bbox_props = dict(facecolor=\"none\", edgecolor=\"none\", lw=0)\n",
"\n",
" if stroke is True:\n",
" patheffects = [PathEffects.withStroke(linewidth=2, foreground=\"w\", alpha=0.8)]\n",
" else:\n",
" patheffects = None\n",
" \n",
" # rotation is 0 for east, rotates CCW\n",
" # rotation is with respect for figure\n",
" ax.annotate(\n",
" txt.format(value=value, axis=axis),\n",
" rotation=fangle_deg + txt_rotation_offset, color=color,\n",
" # rotation_mode=\"default\",\n",
" rotation_mode=\"anchor\",\n",
" bbox=bbox_props, size=12,\n",
" weight=\"bold\", ha=ha, va=va,\n",
" xy=(x0, y0), xycoords='data',\n",
" xytext=(px0, py0), textcoords='offset pixels',\n",
" path_effects=patheffects,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "94c948ca-0b6f-404e-b47c-6dceaba2f052",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"clay_color = \"C3\"\n",
"silt_color = \"C1\"\n",
"sand_color = \"C2\"\n",
"\n",
"fig, ax = plt.subplots(figsize=(10, 6), layout=\"compressed\", facecolor=\"w\", dpi=120)\n",
"\n",
"ax.add_patch(mpl.patches.Polygon(np.stack([boundary_dg, boundary_σg], axis=1), facecolor=\"#F4F4F9\", edgecolor=\"none\", zorder=1, alpha=0.5))\n",
"\n",
"txt_legend1 = \"\"\n",
"txt_legend2 = \"\"\n",
"\n",
"for i, (class_name, (clay, silt)) in enumerate(class_label_pos.items()):\n",
" sand = 100 - (clay+silt)\n",
" x,y = get_geometric_props(clay, silt, sand)\n",
"\n",
" ax.text(\n",
" x, y, class_name.split(\")\")[0], weight=\"bold\", \n",
" color=\"k\", size=14, va=\"center\", ha=\"center\",\n",
" alpha=0.5,\n",
" path_effects=[PathEffects.withStroke(linewidth=2, foreground=\"w\", alpha=.5)],\n",
" )\n",
"\n",
" txt_legend1 += f\"{class_name.split(\" \")[0]}\\n\"\n",
" txt_legend2 += f\"{class_name.title().split(\") \")[1]}\\n\"\n",
" \n",
"ax.text(\n",
" 0.08, 0.9, txt_legend1, va=\"top\", ha=\"right\", \n",
" size=11, weight=\"normal\", alpha=0.8, transform=ax.transAxes,\n",
" path_effects=[PathEffects.withStroke(linewidth=4, foreground=\"w\", alpha=.6)],\n",
")\n",
"\n",
"ax.text(\n",
" 0.085, 0.9, txt_legend2, va=\"top\", ha=\"left\", \n",
" size=11, weight=\"normal\", alpha=0.8, transform=ax.transAxes,\n",
" path_effects=[PathEffects.withStroke(linewidth=4, foreground=\"w\", alpha=.6)],\n",
")\n",
"\n",
"for boundaries, color in [(clay_boundaries, clay_color), (silt_boundaries, silt_color), (sand_boundaries, sand_color)]:\n",
" for dg, σg in boundaries:\n",
" ax.plot(dg, σg, \"-\", lw=.3, alpha=1, color=color)\n",
" \n",
"for name, nodes in class_bounds.items():\n",
" for dg, σg in nodes:\n",
" ax.plot(dg, σg, \"-\", lw=1, alpha=1, color=\"k\")\n",
"\n",
"ax.set_xscale(\"log\", base=10)\n",
"ax.set_xlabel(\"Geometric mean diameter\")\n",
"ax.set_ylabel(\"Geometric\\nstandard deviation\", rotation=0, labelpad=45)\n",
"ax.set_title(\"Soil texture diagram (Shirazi & Boersma (1984))\")\n",
"\n",
"func = get_geometric_props\n",
"\n",
"lbl_min=20\n",
"lbl_max=80\n",
"lbl_step=20\n",
"\n",
"debug = False\n",
"\n",
"annotate(\n",
" ax, label_min=10, label_max=90, label_step=10, position_offset=-15, ha=\"left\", va=\"center\", \n",
" position_rotation_offset=0, txt_rotation_offset=180, color=clay_color, plot_marker=True, axis=\"Clay\", \n",
" txt=\"{value:1.0f}%\", coord_conversion=func, marker_angle=None, align_with=\"grid\", debug=debug, stroke=True,\n",
")\n",
"annotate(\n",
" ax, label_min=lbl_min, label_max=lbl_max, label_step=lbl_step, position_offset=18, ha=\"right\", va=\"center\", \n",
" position_rotation_offset=0, txt_rotation_offset=180, color=silt_color,plot_marker=True, axis=\"Silt\", \n",
" txt=\"{value:1.0f}%\", coord_conversion=func, marker_angle=None, align_with=\"grid\", debug=debug, stroke=True,\n",
")\n",
"annotate(\n",
" ax, label_min=lbl_min, label_max=lbl_max, label_step=lbl_step, position_offset=18, ha=\"left\", va=\"center\", \n",
" position_rotation_offset=0, txt_rotation_offset=0, color=sand_color,plot_marker=True, axis=\"Sand\", \n",
" txt=\"{value:1.0f}%\", coord_conversion=func, marker_angle=None, align_with=\"grid\", debug=debug, stroke=True,\n",
")\n",
"\n",
"annotate(\n",
" ax, label_min=25, label_max=75, label_step=49, position_offset=15, ha=\"center\", va=\"center\", position_rotation_offset=-90, \n",
" txt_rotation_offset=180, color=clay_color, plot_marker=False, axis=\"Clay\", txt=r\"$\\leftarrow${axis}\", \n",
" coord_conversion=func, align_with=\"boundary\", debug=debug,\n",
")\n",
"annotate(\n",
" ax, label_min=50, label_max=50, label_step=50, position_offset=55, ha=\"center\", va=\"center\", position_rotation_offset=-90, \n",
" txt_rotation_offset=0, color=silt_color, plot_marker=False, axis=\"Silt\", txt=r\"{axis}$\\rightarrow$\", \n",
" coord_conversion=func, align_with=\"boundary\", debug=debug,\n",
")\n",
"annotate(\n",
" ax, label_min=50, label_max=50, label_step=50, position_offset=55, ha=\"center\", va=\"center\", position_rotation_offset=-90, \n",
" txt_rotation_offset=0, color=sand_color, plot_marker=False, axis=\"Sand\", txt=r\"{axis}$\\rightarrow$\", \n",
" coord_conversion=func, align_with=\"boundary\", debug=debug,\n",
")\n",
"\n",
"ax.grid(True, linestyle=\"--\", dashes=(10,10), color=\"k\", lw=0.5, alpha=.2, zorder=199)\n",
"ax.set_ylim(-3, 52)\n",
"ax.set_xlim(0.1, 450)\n",
"ax.xaxis.set_major_locator(mpl.ticker.LogLocator(base=10, subs=(1, 2, 5)))\n",
"ax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{round(x, 3)} μm\"))\n",
"ax.set_yticks(range(0, 55, 5))\n",
"ax.set_yticklabels(map(lambda x: f\"{x:1.0f} μm\", ax.get_yticks()))\n",
"ax.set_facecolor(\"w\")\n",
"\n",
"fig.savefig(\"soil_texture_geometric_grid_aligned.png\", dpi=120, bbox_inches=0, pad_inches=0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "39e3835e-ae02-41ca-bb32-43267153558d",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "f5f23f39-876a-4734-b21e-f2abb290a506",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment