Exercise 7.3
Density-Neutron Crossplot
The density-neutron crossplot is how a petrophysicist reads lithology and porosity at a glance. Density porosity (PHID) goes on the x-axis, neutron porosity (NPHI) on the y-axis, and every depth sample becomes a dot coloured by how shaly it is (VSHALE). Clean sand plots near the 1:1 line where PHID ≈ NPHI; shale pulls points up and to the left (high NPHI, low PHID), and gas pulls them the other way (crossover).
You're looking at the OD-003 log over OML 58. Write make_crossplot(logs) that:
- Adds two columns to
logs(modify the frame in place, or work on
a copy you return (the grader reads logs after you call it):
PHID = (2.65 - RHOB) / 1.65, clipped to[0, 0.45]
(matrix density 2.65 g/cc, fluid 1.0 g/cc).
VSHALEfrom the linear gamma-ray index, using
gr_clean = logs["GR"].quantile(0.05) and gr_shale = logs["GR"].quantile(0.95): IGR = (GR - gr_clean) / (gr_shale - gr_clean), then clip to [0, 1].
- Builds and returns a matplotlib
Figurewith a scatter:
x=PHID, y=NPHI, c=VSHALE (a colourmap), plus a 1:1 reference line (ax.plot), an x-label, a y-label, and a title.
Return the Figure object. The grader inspects the axes; it never diffs pixels, so colours and styling are up to you as long as the artists are there.
Stuck? Reveal hints one at a time — they progress from nudge to near-solution.
visibilityReveal reference solutionexpand_more
Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(7)
depth = np.arange(7000.0, 7850.0, 0.5)
_cond = [depth < 7200, depth < 7400, depth < 7600, depth < 7700, depth >= 7700]
_gr = np.select(_cond, [105.0, 40.0, 32.0, 110.0, 45.0])
_rt = np.select(_cond, [4.0, 2.2, 35.0, 4.0, 2.5])
_rhob = np.select(_cond, [2.54, 2.31, 2.26, 2.55, 2.33])
_nphi = np.select(_cond, [0.33, 0.21, 0.17, 0.34, 0.22])
logs = pd.DataFrame({
"DEPTH": depth,
"GR": _gr + np.random.normal(0, 4.0, len(depth)),
"RT": _rt * np.exp(np.random.normal(0, 0.10, len(depth))),
"RHOB": _rhob + np.random.normal(0, 0.02, len(depth)),
"NPHI": np.clip(_nphi + np.random.normal(0, 0.015, len(depth)), 0.0, 1.0),
})
STEP_FT = 0.5
def make_crossplot(logs):
"""Add PHID + VSHALE to `logs`, then return a density-neutron crossplot Figure."""
# 1. Density porosity: matrix 2.65 g/cc, fluid 1.0 g/cc -> denom 1.65.
logs["PHID"] = np.clip((2.65 - logs["RHOB"]) / 1.65, 0.0, 0.45)
# 2. Linear Vshale from GR endpoints (5th / 95th percentile).
gr_clean = logs["GR"].quantile(0.05)
gr_shale = logs["GR"].quantile(0.95)
igr = (logs["GR"] - gr_clean) / (gr_shale - gr_clean)
logs["VSHALE"] = np.clip(igr, 0.0, 1.0)
# 3. Crossplot: PHID (x) vs NPHI (y), coloured by VSHALE, with a 1:1 line.
fig, ax = plt.subplots(figsize=(6, 6))
sc = ax.scatter(logs["PHID"], logs["NPHI"], c=logs["VSHALE"],
cmap="viridis", s=12, alpha=0.8)
ax.plot([0.0, 0.45], [0.0, 0.45], "k--", lw=1.2, label="1:1 (clean sand)")
fig.colorbar(sc, ax=ax, label="VSHALE (v/v)")
ax.set_xlabel("Density porosity PHID (v/v)")
ax.set_ylabel("Neutron porosity NPHI (v/v)")
ax.set_title("OD-003 density-neutron crossplot, OML 58")
ax.legend(loc="upper left")
ax.grid(True, alpha=0.25)
return fig
fig = make_crossplot(logs)
print("PHID range:", round(logs["PHID"].min(), 3), "->", round(logs["PHID"].max(), 3))
print("VSHALE range:", round(logs["VSHALE"].min(), 3), "->", round(logs["VSHALE"].max(), 3))
plt.show()
lockCopying code is a Full Access feature.