|
| 1 | +""" pyplots.ai |
| 2 | +column-stratigraphic: Stratigraphic Column with Lithology Patterns |
| 3 | +Library: letsplot 4.9.0 | Python 3.14.3 |
| 4 | +Quality: 88/100 | Created: 2026-03-15 |
| 5 | +""" |
| 6 | +# ruff: noqa: F405 |
| 7 | + |
| 8 | +import os |
| 9 | +import shutil |
| 10 | + |
| 11 | +import pandas as pd |
| 12 | +from lets_plot import * # noqa: F403 |
| 13 | + |
| 14 | + |
| 15 | +LetsPlot.setup_html() |
| 16 | + |
| 17 | +# Data - synthetic sedimentary section with 10 layers (depth increasing downward) |
| 18 | +layers = pd.DataFrame( |
| 19 | + { |
| 20 | + "top": [0, 15, 35, 55, 80, 110, 135, 165, 195, 230], |
| 21 | + "bottom": [15, 35, 55, 80, 110, 135, 165, 195, 230, 260], |
| 22 | + "lithology": [ |
| 23 | + "Sandstone", |
| 24 | + "Shale", |
| 25 | + "Limestone", |
| 26 | + "Siltstone", |
| 27 | + "Sandstone", |
| 28 | + "Conglomerate", |
| 29 | + "Shale", |
| 30 | + "Limestone", |
| 31 | + "Siltstone", |
| 32 | + "Sandstone", |
| 33 | + ], |
| 34 | + "formation": [ |
| 35 | + "Basal Sand Fm", |
| 36 | + "Dark Shale Mbr", |
| 37 | + "Reef Limestone Fm", |
| 38 | + "Gray Silt Mbr", |
| 39 | + "Channel Sand Fm", |
| 40 | + "Gravel Bed Fm", |
| 41 | + "Marine Shale Fm", |
| 42 | + "Platform Carb Fm", |
| 43 | + "Tidal Flat Mbr", |
| 44 | + "Upper Sand Fm", |
| 45 | + ], |
| 46 | + "age": [ |
| 47 | + "Triassic", |
| 48 | + "Triassic", |
| 49 | + "Jurassic", |
| 50 | + "Jurassic", |
| 51 | + "Jurassic", |
| 52 | + "Cretaceous", |
| 53 | + "Cretaceous", |
| 54 | + "Cretaceous", |
| 55 | + "Paleogene", |
| 56 | + "Paleogene", |
| 57 | + ], |
| 58 | + } |
| 59 | +) |
| 60 | + |
| 61 | +layers["thickness"] = layers["bottom"] - layers["top"] |
| 62 | +layers["xmin"] = 0.0 |
| 63 | +layers["xmax"] = 1.0 |
| 64 | + |
| 65 | +# Lithology color palette (refined earth tones for geological context) |
| 66 | +lithology_colors = { |
| 67 | + "Sandstone": "#F2D16B", |
| 68 | + "Shale": "#8C8C8C", |
| 69 | + "Limestone": "#6AAED6", |
| 70 | + "Siltstone": "#C08C5A", |
| 71 | + "Conglomerate": "#D06A62", |
| 72 | +} |
| 73 | +lithology_order = ["Sandstone", "Shale", "Limestone", "Siltstone", "Conglomerate"] |
| 74 | + |
| 75 | +# Generate lithology pattern overlay data |
| 76 | +pattern_segments = [] # For line-based patterns (shale dashes, limestone bricks) |
| 77 | +pattern_points = [] # For dot-based patterns (sandstone stipple, conglomerate circles) |
| 78 | + |
| 79 | +for _, row in layers.iterrows(): |
| 80 | + top, bottom, lith = row["top"], row["bottom"], row["lithology"] |
| 81 | + thickness = bottom - top |
| 82 | + margin = 0.03 |
| 83 | + |
| 84 | + if lith == "Shale": |
| 85 | + # Horizontal dashes - denser pattern |
| 86 | + spacing = 2.5 |
| 87 | + n_lines = max(1, int(thickness / spacing)) |
| 88 | + for i in range(n_lines): |
| 89 | + y = top + (i + 0.5) * thickness / n_lines |
| 90 | + for x_start in [0.04, 0.20, 0.36, 0.52, 0.68, 0.84]: |
| 91 | + pattern_segments.append({"x": x_start, "y": y, "xend": x_start + 0.12, "yend": y}) |
| 92 | + |
| 93 | + elif lith == "Limestone": |
| 94 | + # Brick pattern: horizontal lines + offset vertical lines - denser |
| 95 | + spacing = 3.5 |
| 96 | + n_rows = max(1, int(thickness / spacing)) |
| 97 | + for i in range(n_rows + 1): |
| 98 | + y = top + margin + i * (thickness - 2 * margin) / max(n_rows, 1) |
| 99 | + if top + margin <= y <= bottom - margin: |
| 100 | + pattern_segments.append({"x": 0.02, "y": y, "xend": 0.98, "yend": y}) |
| 101 | + for i in range(n_rows): |
| 102 | + y_top = top + margin + i * (thickness - 2 * margin) / max(n_rows, 1) |
| 103 | + y_bot = top + margin + (i + 1) * (thickness - 2 * margin) / max(n_rows, 1) |
| 104 | + offset = 0.2 if i % 2 == 0 else 0.0 |
| 105 | + for vx in [0.2 + offset, 0.4 + offset, 0.6 + offset, 0.8 + offset]: |
| 106 | + if 0.02 < vx < 0.98: |
| 107 | + pattern_segments.append({"x": vx, "y": y_top, "xend": vx, "yend": y_bot}) |
| 108 | + |
| 109 | + elif lith == "Sandstone": |
| 110 | + # Stipple dots - denser grid |
| 111 | + spacing_y = 2.8 |
| 112 | + spacing_x = 0.08 |
| 113 | + n_rows = max(1, int(thickness / spacing_y)) |
| 114 | + for i in range(n_rows): |
| 115 | + y = top + (i + 0.5) * thickness / n_rows |
| 116 | + offset = 0.04 if i % 2 == 0 else 0.0 |
| 117 | + x = 0.06 + offset |
| 118 | + while x < 0.96: |
| 119 | + pattern_points.append({"x": x, "y": y, "shape": "dot"}) |
| 120 | + x += spacing_x |
| 121 | + |
| 122 | + elif lith == "Siltstone": |
| 123 | + # Short tilted dashes - denser |
| 124 | + spacing_y = 3.5 |
| 125 | + n_rows = max(1, int(thickness / spacing_y)) |
| 126 | + for i in range(n_rows): |
| 127 | + y = top + (i + 0.5) * thickness / n_rows |
| 128 | + offset = 0.06 if i % 2 == 0 else 0.0 |
| 129 | + for x in [0.08 + offset, 0.22 + offset, 0.36 + offset, 0.50 + offset, 0.64 + offset, 0.78 + offset]: |
| 130 | + if x < 0.95: |
| 131 | + pattern_segments.append({"x": x, "y": y - 0.8, "xend": x + 0.06, "yend": y + 0.8}) |
| 132 | + |
| 133 | + elif lith == "Conglomerate": |
| 134 | + # Circles (larger dots) - denser |
| 135 | + spacing_y = 4.5 |
| 136 | + n_rows = max(1, int(thickness / spacing_y)) |
| 137 | + for i in range(n_rows): |
| 138 | + y = top + (i + 0.5) * thickness / n_rows |
| 139 | + offset = 0.08 if i % 2 == 0 else 0.0 |
| 140 | + for x in [0.12 + offset, 0.32 + offset, 0.52 + offset, 0.72 + offset]: |
| 141 | + if x < 0.95: |
| 142 | + pattern_points.append({"x": x, "y": y, "shape": "circle"}) |
| 143 | + |
| 144 | +pattern_seg_df = pd.DataFrame(pattern_segments) if pattern_segments else None |
| 145 | +pattern_dot_df = pd.DataFrame([p for p in pattern_points if p["shape"] == "dot"]) if pattern_points else None |
| 146 | +pattern_circle_df = pd.DataFrame([p for p in pattern_points if p["shape"] == "circle"]) if pattern_points else None |
| 147 | + |
| 148 | +# Formation labels (right side) |
| 149 | +form_labels = [] |
| 150 | +for _, row in layers.iterrows(): |
| 151 | + mid_depth = (row["top"] + row["bottom"]) / 2 |
| 152 | + form_labels.append({"x": 1.08, "y": mid_depth, "label": row["formation"]}) |
| 153 | +form_df = pd.DataFrame(form_labels) |
| 154 | + |
| 155 | +# Age labels (left side) with bracket indicators |
| 156 | +age_spans = {"Triassic": (0, 35), "Jurassic": (35, 110), "Cretaceous": (110, 195), "Paleogene": (195, 260)} |
| 157 | +age_labels = [] |
| 158 | +age_brackets = [] |
| 159 | +for age_name, (age_top, age_bottom) in age_spans.items(): |
| 160 | + age_labels.append({"x": -0.18, "y": (age_top + age_bottom) / 2, "label": age_name}) |
| 161 | + # Bracket lines on left |
| 162 | + age_brackets.append({"x": -0.06, "y": age_top + 1, "xend": -0.06, "yend": age_bottom - 1}) |
| 163 | + age_brackets.append({"x": -0.06, "y": age_top + 1, "xend": -0.03, "yend": age_top + 1}) |
| 164 | + age_brackets.append({"x": -0.06, "y": age_bottom - 1, "xend": -0.03, "yend": age_bottom - 1}) |
| 165 | +age_df = pd.DataFrame(age_labels) |
| 166 | +bracket_df = pd.DataFrame(age_brackets) |
| 167 | + |
| 168 | +# Unconformity wavy line at Jurassic/Cretaceous boundary (110m) |
| 169 | +wavy_x = [] |
| 170 | +wavy_y = [] |
| 171 | +n_waves = 20 |
| 172 | +for i in range(n_waves + 1): |
| 173 | + xi = i / n_waves |
| 174 | + yi = 110 + 1.5 * (1 if (i % 2 == 0) else -1) |
| 175 | + wavy_x.append(xi) |
| 176 | + wavy_y.append(yi) |
| 177 | +wavy_df = pd.DataFrame({"x": wavy_x, "y": wavy_y}) |
| 178 | + |
| 179 | +# Plot assembly |
| 180 | +plot = ( |
| 181 | + ggplot() |
| 182 | + # Layer rectangles with interactive tooltips |
| 183 | + + geom_rect( |
| 184 | + aes(xmin="xmin", xmax="xmax", ymin="top", ymax="bottom", fill="lithology"), |
| 185 | + data=layers, |
| 186 | + color="#2C2C2C", |
| 187 | + size=1.0, |
| 188 | + alpha=0.8, |
| 189 | + tooltips=layer_tooltips() |
| 190 | + .format("@top", ".0f") |
| 191 | + .format("@bottom", ".0f") |
| 192 | + .format("@thickness", ".0f") |
| 193 | + .title("@formation") |
| 194 | + .line("@lithology | @age") |
| 195 | + .line("Depth: @top\u2013@bottom m") |
| 196 | + .line("Thickness: @thickness m"), |
| 197 | + ) |
| 198 | +) |
| 199 | + |
| 200 | +# Add pattern overlays |
| 201 | +if pattern_seg_df is not None and len(pattern_seg_df) > 0: |
| 202 | + plot = plot + geom_segment( |
| 203 | + aes(x="x", y="y", xend="xend", yend="yend"), |
| 204 | + data=pattern_seg_df, |
| 205 | + color="#2A2A2A", |
| 206 | + size=0.6, |
| 207 | + alpha=0.75, |
| 208 | + show_legend=False, |
| 209 | + ) |
| 210 | + |
| 211 | +if pattern_dot_df is not None and len(pattern_dot_df) > 0: |
| 212 | + plot = plot + geom_point( |
| 213 | + aes(x="x", y="y"), data=pattern_dot_df, color="#4A3A10", size=1.8, alpha=0.7, shape=16, show_legend=False |
| 214 | + ) |
| 215 | + |
| 216 | +if pattern_circle_df is not None and len(pattern_circle_df) > 0: |
| 217 | + plot = plot + geom_point( |
| 218 | + aes(x="x", y="y"), data=pattern_circle_df, color="#5A2A2A", size=5.0, alpha=0.7, shape=1, show_legend=False |
| 219 | + ) |
| 220 | + |
| 221 | +# Unconformity wavy line at 110m with label |
| 222 | +plot = plot + geom_line(aes(x="x", y="y"), data=wavy_df, color="#C0392B", size=2.0, show_legend=False) |
| 223 | +unconformity_label = pd.DataFrame({"x": [1.08], "y": [110], "label": ["Unconformity"]}) |
| 224 | +plot = plot + geom_text( |
| 225 | + aes(x="x", y="y", label="label"), |
| 226 | + data=unconformity_label, |
| 227 | + color="#C0392B", |
| 228 | + size=13, |
| 229 | + fontface="bold", |
| 230 | + hjust=0, |
| 231 | + show_legend=False, |
| 232 | +) |
| 233 | + |
| 234 | +# Age boundary dashed lines (non-unconformity) |
| 235 | +non_unconformity_boundaries = pd.DataFrame( |
| 236 | + {"x": [-0.02, -0.02], "y": [35, 195], "xend": [1.02, 1.02], "yend": [35, 195]} |
| 237 | +) |
| 238 | +plot = plot + geom_segment( |
| 239 | + aes(x="x", y="y", xend="xend", yend="yend"), |
| 240 | + data=non_unconformity_boundaries, |
| 241 | + linetype="dashed", |
| 242 | + color="#666666", |
| 243 | + size=0.6, |
| 244 | + show_legend=False, |
| 245 | +) |
| 246 | + |
| 247 | +# Age brackets (left side) |
| 248 | +plot = plot + geom_segment( |
| 249 | + aes(x="x", y="y", xend="xend", yend="yend"), data=bracket_df, color="#444444", size=0.6, show_legend=False |
| 250 | +) |
| 251 | + |
| 252 | +# Formation labels (right side) |
| 253 | +plot = plot + geom_text(aes(x="x", y="y", label="label"), data=form_df, size=15, color="#2C2C2C", hjust=0) |
| 254 | + |
| 255 | +# Age labels (left side) |
| 256 | +plot = plot + geom_text(aes(x="x", y="y", label="label"), data=age_df, size=15, color="#2C2C2C", fontface="italic") |
| 257 | + |
| 258 | +# Scales and theme |
| 259 | +plot = ( |
| 260 | + plot |
| 261 | + + scale_fill_manual(values=lithology_colors, name="Lithology", limits=lithology_order) |
| 262 | + + scale_y_reverse() |
| 263 | + + labs( |
| 264 | + title="column-stratigraphic \u00b7 letsplot \u00b7 pyplots.ai", |
| 265 | + subtitle="Synthetic Mesozoic\u2013Cenozoic sedimentary section \u00b7 Triassic to Paleogene", |
| 266 | + y="Depth (m)", |
| 267 | + x="", |
| 268 | + ) |
| 269 | + + scale_x_continuous(limits=[-0.28, 1.52]) |
| 270 | + + flavor_high_contrast_light() |
| 271 | + + theme( |
| 272 | + plot_title=element_text(size=26, face="bold", color="#1A1A1A", margin=[0, 0, 12, 0]), |
| 273 | + plot_subtitle=element_text(size=16, color="#666666", face="italic"), |
| 274 | + axis_title_y=element_text(size=20, color="#333333", margin=[0, 8, 0, 0]), |
| 275 | + axis_title_x=element_blank(), |
| 276 | + axis_text_y=element_text(size=16, color="#444444"), |
| 277 | + axis_text_x=element_blank(), |
| 278 | + axis_ticks_x=element_blank(), |
| 279 | + axis_line_y=element_line(size=0.8, color="#333333"), |
| 280 | + legend_title=element_text(size=16, face="bold"), |
| 281 | + legend_text=element_text(size=14), |
| 282 | + legend_position="bottom", |
| 283 | + panel_grid_major_x=element_blank(), |
| 284 | + panel_grid_minor_x=element_blank(), |
| 285 | + panel_grid_major_y=element_line(size=0.2, color="#D8D8D8"), |
| 286 | + panel_grid_minor_y=element_blank(), |
| 287 | + plot_background=element_rect(color="white", fill="white"), |
| 288 | + plot_margin=[16, 16, 16, 16], |
| 289 | + ) |
| 290 | + + ggsize(1600, 900) |
| 291 | +) |
| 292 | + |
| 293 | +# Save |
| 294 | +ggsave(plot, "plot.png", scale=3) |
| 295 | +ggsave(plot, "plot.html") |
| 296 | + |
| 297 | +# Move files from lets-plot-images to current directory |
| 298 | +if os.path.exists("lets-plot-images/plot.png"): |
| 299 | + shutil.move("lets-plot-images/plot.png", "plot.png") |
| 300 | +if os.path.exists("lets-plot-images/plot.html"): |
| 301 | + shutil.move("lets-plot-images/plot.html", "plot.html") |
| 302 | +if os.path.exists("lets-plot-images") and not os.listdir("lets-plot-images"): |
| 303 | + os.rmdir("lets-plot-images") |
0 commit comments