Lucas e' un database europeo di analisi di alcuni parametri di suoli collegati agli spettri MIR e VNIR degli stessi campioni georiferiti
Le analisi di laboratorio sono per i parametri sottostanti (a fianco il numero di campioni analizzati...non tuti campioni del database ha lo stessa caratterizzazione chimica)
Il sito originale del progetto Luca e' ma risulta piu' comodo https://docs.soilspectroscopy.org/datasets/LUCAS.html dove non e' richiesta l'autorizzazione
| Scheletro (frammenti grossolani) | % peso | ~26.300 |
| Argilla | % peso | ~26.300 |
| Limo | % peso | ~26.300 |
| Sabbia | % peso | ~26.300 |
| pH in H2O | indice | ~43.900 |
| pH in CaCl2 | indice | ~43.900 |
| Carbonio organico | % peso | ~43.700 |
| CaCO3 | % peso | ~43.200 |
| Azoto totale | % peso | ~43.900 |
| Fosforo estraibile | mg/kg | ~43.200 |
| Potassio estraibile | cmol/kg | ~43.900 |
| CEC | cmol/kg | ~22.000 |
| Conducibilità elettrica | dS/m | ~21.900 |
#!/usr/bin/env python3
"""
Estrae da LUCAS (OSSL v1.3, https://docs.soilspectroscopy.org/datasets/LUCAS.html)
un CSV con: latitudine, longitudine, CaCO3 (w.pct) e spettro VNIR (riflettanza 400-2500 nm).
Uso:
pip install pandas pyarrow requests
python lucas_caco3_vnir.py # -> lucas_caco3_vnir.csv
python lucas_caco3_vnir.py -o out.csv --bbox 6.6 35.5 18.6 47.1 # solo Italia (lon_min lat_min lon_max lat_max)
python lucas_caco3_vnir.py --country IT # filtro per paese (se la colonna e' popolata)
"""
import argparse
import re
import sys
from pathlib import Path
import pandas as pd
import pyarrow.parquet as pq
import requests
BASE = "https://storage.googleapis.com/soilspec4gg-public/datasets/LUCAS/"
FILES = {
"site": "ossl_soilsite_v1.3.parquet",
"lab": "ossl_soillab_v1.3.parquet",
"visnir": "ossl_visnir_v1.3.parquet",
}
ID = "id.layer_uuid_txt"
LON = "longitude.point_wgs84_dd"
LAT = "latitude.point_wgs84_dd"
COUNTRY = "loc.country_src_txt"
CACO3 = "caco3_iso.10693_w.pct"
def download(name: str, cache: Path) -> Path:
cache.mkdir(parents=True, exist_ok=True)
dest = cache / name
if dest.exists() and dest.stat().st_size > 0:
return dest
url = BASE + name
print(f"Scarico {url}", file=sys.stderr)
tmp = dest.with_suffix(dest.suffix + ".part")
with requests.get(url, stream=True, timeout=60) as r:
r.raise_for_status()
total = int(r.headers.get("content-length", 0))
done = 0
with open(tmp, "wb") as f:
for chunk in r.iter_content(chunk_size=1 << 20):
f.write(chunk)
done += len(chunk)
if total:
print(f"\r {done / total:6.1%}", end="", file=sys.stderr)
print(file=sys.stderr)
tmp.rename(dest)
return dest
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("-o", "--output", default="lucas_caco3_vnir.csv")
ap.add_argument("--cache", default="./ossl_cache", help="cartella per i file scaricati")
ap.add_argument("--bbox", nargs=4, type=float, metavar=("LON_MIN", "LAT_MIN", "LON_MAX", "LAT_MAX"))
ap.add_argument("--country", help="codice paese (es. IT), usa la colonna loc.country_src_txt")
args = ap.parse_args()
cache = Path(args.cache)
# --- coordinate (e paese, se presente)
site_path = download(FILES["site"], cache)
site_cols = [c for c in (ID, LON, LAT, COUNTRY) if c in pq.ParquetFile(site_path).schema.names]
site = pd.read_parquet(site_path, columns=site_cols)
# --- CaCO3
lab_path = download(FILES["lab"], cache)
lab_names = pq.ParquetFile(lab_path).schema.names
caco3_col = CACO3 if CACO3 in lab_names else next((c for c in lab_names if c.startswith("caco3")), None)
if caco3_col is None:
sys.exit("Colonna CaCO3 non trovata nel file soillab.")
lab = pd.read_parquet(lab_path, columns=[ID, caco3_col]).rename(columns={caco3_col: "caco3_wpct"})
# --- spettri VNIR (solo colonne di riflettanza)
vn_path = download(FILES["visnir"], cache)
vn_names = pq.ParquetFile(vn_path).schema.names
spec_cols = [c for c in vn_names if c.startswith("scan_visnir") and c.endswith("_ref")]
if not spec_cols:
sys.exit("Nessuna colonna scan_visnir.*_ref trovata nel file visnir.")
# ordina per lunghezza d'onda e rinomina in 'nm'
wl = {c: int(float(re.search(r"scan_visnir\.([\d.]+)_ref", c).group(1))) for c in spec_cols}
spec_cols.sort(key=wl.get)
vn = pd.read_parquet(vn_path, columns=[ID] + spec_cols).rename(columns=wl)
# --- join sull'id del campione
df = site.merge(lab, on=ID, how="inner").merge(vn, on=ID, how="inner")
# --- filtri
if args.country and COUNTRY in df:
df = df[df[COUNTRY].astype(str).str.upper() == args.country.upper()]
if args.bbox:
lon0, lat0, lon1, lat1 = args.bbox
df = df[df[LON].between(lon0, lon1) & df[LAT].between(lat0, lat1)]
df = df.dropna(subset=[LON, LAT, "caco3_wpct"])
# --- output
out = df.rename(columns={LON: "longitude", LAT: "latitude"})
cols = ["latitude", "longitude", "caco3_wpct"] + list(wl.values())
out[cols].to_csv(args.output, index=False)
print(f"Scritte {len(out)} righe, {len(wl)} bande ({min(wl.values())}-{max(wl.values())} nm) -> {args.output}")
if __name__ == "__main__":
main()
Per questa prova lo script sovrastante scarica solo la posizione, il valore di laboratorio di CaCO3 e la firma spettrale VNIRSWIR e salva in un file csv
La correlazione tra i parametri e' invece fatta tramite
#!/usr/bin/env python3
"""
Correlazioni tra le proprieta' di laboratorio del dataset LUCAS (OSSL soillab).
Calcola, con campioni a coppie completi (pairwise), la correlazione di Spearman (consigliata: le
distribuzioni sono molto asimmetriche) e di Pearson, il numero di campioni per coppia e un heatmap.
Uso:
pip install pandas pyarrow requests matplotlib
python lucas_correlations.py # scarica soillab in ./ossl_cache se manca
python lucas_correlations.py --parquet ossl_cache/ossl_soillab_v1.3.parquet --log
Output (prefisso --prefix, default lucas_corr):
<prefix>_spearman.csv, <prefix>_pearson.csv, <prefix>_n.csv, <prefix>_pairs.csv, <prefix>_heatmap.png
"""
import argparse
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
URL = "https://storage.googleapis.com/soilspec4gg-public/datasets/LUCAS/ossl_soillab_v1.3.parquet"
SKEWED = ["oc", "caco3", "n.tot", "p.ext", "k.ext", "cec", "ec", "cf"] # variabili con coda lunga
def get_file(path: Path) -> Path:
if path.exists() and path.stat().st_size > 0:
return path
path.parent.mkdir(parents=True, exist_ok=True)
print(f"Scarico {URL}")
with requests.get(URL, stream=True, timeout=60) as r:
r.raise_for_status()
with open(path, "wb") as f:
for chunk in r.iter_content(1 << 20):
f.write(chunk)
return path
def short(col: str) -> str:
"""'clay.tot_iso.11277_w.pct' -> 'clay.tot'"""
return col.split("_")[0]
def heatmap(corr: pd.DataFrame, path: str, title: str):
n = len(corr)
fig, ax = plt.subplots(figsize=(0.75 * n + 3, 0.75 * n + 2.5))
masked = np.ma.masked_where(np.triu(np.ones_like(corr, dtype=bool), k=1), corr.to_numpy())
im = ax.imshow(masked, cmap="RdBu_r", vmin=-1, vmax=1)
ax.set_xticks(range(n), corr.columns, rotation=45, ha="right")
ax.set_yticks(range(n), corr.index)
for i in range(n):
for j in range(i + 1):
v = corr.iat[i, j]
if np.isfinite(v):
ax.text(j, i, f"{v:.2f}", ha="center", va="center", fontsize=8,
color="white" if abs(v) > 0.6 else "black")
fig.colorbar(im, ax=ax, shrink=0.7, label="correlazione")
ax.set_title(title)
fig.tight_layout()
fig.savefig(path, dpi=150)
plt.close(fig)
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--parquet", default="ossl_cache/ossl_soillab_v1.3.parquet")
ap.add_argument("--prefix", default="lucas_corr")
ap.add_argument("--log", action="store_true", help="Pearson su log1p delle variabili asimmetriche")
ap.add_argument("--min-n", type=int, default=30, help="campioni minimi per calcolare una coppia")
args = ap.parse_args()
df = pd.read_parquet(get_file(Path(args.parquet)))
cols = [c for c in df.columns if not c.startswith(("id.", "__")) and pd.api.types.is_numeric_dtype(df[c])]
df = df[cols].rename(columns=short)
print(f"{len(df)} campioni, {df.shape[1]} proprieta': {', '.join(df.columns)}\n")
avail = df.notna().sum()
print("Campioni disponibili per proprieta':")
print(avail.to_string(), "\n")
print(f"Campioni con TUTTE le proprieta' presenti: {int(df.notna().all(axis=1).sum())}\n")
n_pair = df.notna().astype(int).T @ df.notna().astype(int)
spearman = df.corr(method="spearman", min_periods=args.min_n)
dfp = df.copy()
if args.log:
for c in SKEWED:
if c in dfp:
dfp[c] = np.log1p(dfp[c].clip(lower=0))
pearson = dfp.corr(method="pearson", min_periods=args.min_n)
spearman.to_csv(f"{args.prefix}_spearman.csv")
pearson.to_csv(f"{args.prefix}_pearson.csv")
n_pair.to_csv(f"{args.prefix}_n.csv")
rows = []
for i, a in enumerate(df.columns):
for b in df.columns[i + 1:]:
rows.append({"var1": a, "var2": b, "spearman": spearman.at[a, b], "pearson": pearson.at[a, b],
"n": int(n_pair.at[a, b])})
pairs = pd.DataFrame(rows)
pairs["abs_spearman"] = pairs["spearman"].abs()
pairs = pairs.sort_values("abs_spearman", ascending=False).drop(columns="abs_spearman")
pairs.to_csv(f"{args.prefix}_pairs.csv", index=False)
pd.set_option("display.width", 120)
print("Coppie ordinate per |Spearman|:")
print(pairs.round(2).to_string(index=False))
heatmap(spearman, f"{args.prefix}_heatmap.png", "LUCAS - correlazione di Spearman (coppie complete)")
print(f"\nSalvati: {args.prefix}_heatmap.png, _spearman.csv, _pearson.csv, _n.csv, _pairs.csv")
if __name__ == "__main__":
main()
Questo e' invece il file che effettua la rete MLP e Random Forest
#!/usr/bin/env python3
"""
Predizione del CaCO3 da spettri VNIR (LUCAS) con Random Forest e rete neurale (MLP).
Split: 70% train / 10% validation / 20% test (stratificato sui quantili di CaCO3).
- validation: scelta degli iperparametri (RF: griglia piccola; MLP: architettura + early stopping)
- test: usato una sola volta per le metriche finali.
Input : CSV con colonne latitude, longitude, caco3_wpct e colonne numeriche = lunghezze d'onda (nm)
Output: <prefix>_performance.png, <prefix>_metrics.csv, <prefix>_predictions.csv
Uso:
pip install numpy pandas scikit-learn scipy matplotlib
python caco3_rf_nn.py lucas_caco3_vnir.csv
python caco3_rf_nn.py lucas_caco3_vnir.csv --band-step 4 --preproc sg1 --log
"""
import argparse
import copy
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.signal import savgol_filter
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import ParameterGrid, train_test_split
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
TARGET = "caco3_wpct"
COLORS = {"train": "tab:blue", "val": "tab:orange", "test": "tab:red"}
# ----------------------------------------------------------------- dati
def is_number(s):
try:
float(s)
return True
except ValueError:
return False
def load(path, band_step):
df = pd.read_csv(path)
spec = [c for c in df.columns if is_number(c)][::band_step]
if not spec:
raise SystemExit("Nessuna colonna spettrale numerica trovata.")
df = df.dropna(subset=[TARGET] + spec).reset_index(drop=True)
X = df[spec].to_numpy(float)
y = df[TARGET].to_numpy(float)
wl = np.array([float(c) for c in spec])
return X, y, wl
def preprocess(X, mode):
if mode == "snv":
return (X - X.mean(1, keepdims=True)) / (X.std(1, keepdims=True) + 1e-12)
if mode == "sg1": # derivata prima Savitzky-Golay
return savgol_filter(X, 11, 2, deriv=1, axis=1)
return X
def split_70_10_20(n, y, seed):
"""Indici train/val/test = 70/10/20, stratificati sui decili di y."""
idx = np.arange(n)
strata = pd.qcut(pd.Series(y).rank(method="first"), 10, labels=False).to_numpy()
rest, test = train_test_split(idx, test_size=0.20, random_state=seed, stratify=strata)
train, val = train_test_split(rest, test_size=0.125, random_state=seed, stratify=strata[rest]) # 0.125*0.8 = 10%
return train, val, test
def metrics(y, p):
rmse = float(np.sqrt(mean_squared_error(y, p)))
iqr = float(np.subtract(*np.percentile(y, [75, 25])))
return {
"n": len(y),
"R2": r2_score(y, p),
"RMSE": rmse,
"MAE": mean_absolute_error(y, p),
"RPIQ": iqr / rmse if rmse > 0 else np.nan,
}
# ----------------------------------------------------------------- modelli
def fit_rf(Xtr, ytr, Xva, yva, inv, n_estimators, seed):
grid = ParameterGrid({"max_features": ["sqrt", 0.2, 0.5], "min_samples_leaf": [1, 3]})
best, best_rmse = None, np.inf
for params in grid:
rf = RandomForestRegressor(n_estimators=n_estimators, n_jobs=-1, random_state=seed, **params)
rf.fit(Xtr, ytr)
rmse = np.sqrt(mean_squared_error(inv(yva), inv(rf.predict(Xva))))
print(f" RF {params} -> val RMSE {rmse:.3f}")
if rmse < best_rmse:
best, best_rmse, best_params = rf, rmse, params
print(f" RF migliore: {best_params}")
return best
def fit_mlp(Xtr, ytr, Xva, yva, seed, epochs=300, patience=30):
"""Prova alcune architetture; early stopping sul RMSE di validation (spazio standardizzato)."""
best = {"rmse": np.inf}
for hidden in [(128,), (256, 64), (512, 128, 32)]:
m = MLPRegressor(hidden_layer_sizes=hidden, alpha=1e-3, solver="adam", learning_rate_init=1e-3,
batch_size=min(64, len(Xtr)), random_state=seed)
run_best, wait, tr_hist, va_hist = np.inf, 0, [], []
snapshot = None
for _ in range(epochs):
m.partial_fit(Xtr, ytr)
tr = np.sqrt(mean_squared_error(ytr, m.predict(Xtr)))
va = np.sqrt(mean_squared_error(yva, m.predict(Xva)))
tr_hist.append(tr)
va_hist.append(va)
if va < run_best - 1e-5:
run_best, wait, snapshot = va, 0, copy.deepcopy(m)
else:
wait += 1
if wait >= patience:
break
print(f" MLP {hidden} -> val RMSE (scalato) {run_best:.3f} in {len(va_hist)} epoche")
if run_best < best["rmse"]:
best = {"rmse": run_best, "model": snapshot, "hidden": hidden, "tr": tr_hist, "va": va_hist}
print(f" MLP migliore: {best['hidden']}")
return best
# ----------------------------------------------------------------- grafici
def scatter(ax, name, preds, ys, mets):
lim = max(max(ys[s].max(), preds[s].max()) for s in ys) * 1.05
ax.plot([0, lim], [0, lim], "k--", lw=1)
for s in ("train", "val", "test"):
ax.scatter(ys[s], preds[s], s=8 if s == "train" else 14, alpha=0.35 if s == "train" else 0.7,
color=COLORS[s], label=f"{s} (R²={mets[s]['R2']:.2f})")
t = mets["test"]
ax.set_title(f"{name} | test: R²={t['R2']:.2f}, RMSE={t['RMSE']:.2f}, RPIQ={t['RPIQ']:.2f}")
ax.set_xlabel("CaCO3 misurato (w%)")
ax.set_ylabel("CaCO3 predetto (w%)")
ax.set_xlim(0, lim)
ax.set_ylim(0, lim)
ax.legend(loc="upper left", fontsize=8)
ax.set_aspect("equal", adjustable="box")
def make_figure(path, preds, ys, mets, rf, wl, mlp):
fig, axes = plt.subplots(2, 2, figsize=(13, 11))
scatter(axes[0, 0], "Random Forest", preds["RF"], ys, mets["RF"])
scatter(axes[0, 1], "Rete neurale (MLP)", preds["MLP"], ys, mets["MLP"])
ax = axes[1, 0]
ax.plot(wl, rf.feature_importances_, color="tab:green")
ax.set_xlabel("Lunghezza d'onda (nm)")
ax.set_ylabel("Importanza (RF)")
ax.set_title("Importanza delle bande VNIR")
ax = axes[1, 1]
ax.plot(mlp["tr"], label="train")
ax.plot(mlp["va"], label="validation")
ax.axvline(int(np.argmin(mlp["va"])), color="gray", ls=":", label="best epoch")
ax.set_xlabel("Epoca")
ax.set_ylabel("RMSE (spazio standardizzato)")
ax.set_title(f"Curva di apprendimento MLP {mlp['hidden']}")
ax.legend()
fig.tight_layout()
fig.savefig(path, dpi=150)
plt.close(fig)
# ----------------------------------------------------------------- main
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("csv")
ap.add_argument("--prefix", default="caco3")
ap.add_argument("--band-step", type=int, default=2, help="usa 1 banda ogni N (default 2 = ogni 4 nm)")
ap.add_argument("--preproc", choices=["none", "snv", "sg1"], default="none")
ap.add_argument("--log", action="store_true", help="modella log1p(CaCO3) (utile: distribuzione molto asimmetrica)")
ap.add_argument("--n-estimators", type=int, default=300)
ap.add_argument("--seed", type=int, default=42)
args = ap.parse_args()
X, y, wl = load(args.csv, args.band_step)
X = preprocess(X, args.preproc)
print(f"{len(y)} campioni, {X.shape[1]} bande ({wl.min():.0f}-{wl.max():.0f} nm)")
tr, va, te = split_70_10_20(len(y), y, args.seed)
print(f"train {len(tr)} ({len(tr)/len(y):.0%}) | val {len(va)} ({len(va)/len(y):.0%}) | test {len(te)} ({len(te)/len(y):.0%})")
sets = {"train": tr, "val": va, "test": te}
fwd = np.log1p if args.log else (lambda v: v)
inv = (lambda v: np.clip(np.expm1(v), 0, None)) if args.log else (lambda v: np.clip(v, 0, None))
ys = {s: y[i] for s, i in sets.items()} # unita' originali
yt = {s: fwd(y[i]) for s, i in sets.items()} # unita' del modello
# --- Random Forest (nessuno scaling necessario)
print("Random Forest...")
rf = fit_rf(X[tr], yt["train"], X[va], yt["val"], inv, args.n_estimators, args.seed)
# --- MLP (X e y standardizzati sul solo train)
print("MLP...")
xs, ysc = StandardScaler().fit(X[tr]), StandardScaler().fit(yt["train"].reshape(-1, 1))
Xs = {s: xs.transform(X[i]) for s, i in sets.items()}
Ys = {s: ysc.transform(yt[s].reshape(-1, 1)).ravel() for s in sets}
mlp = fit_mlp(Xs["train"], Ys["train"], Xs["val"], Ys["val"], args.seed)
# --- predizioni in unita' originali
preds = {"RF": {}, "MLP": {}}
for s, i in sets.items():
preds["RF"][s] = inv(rf.predict(X[i]))
preds["MLP"][s] = inv(ysc.inverse_transform(mlp["model"].predict(Xs[s]).reshape(-1, 1)).ravel())
mets = {m: {s: metrics(ys[s], preds[m][s]) for s in sets} for m in preds}
rows = [{"model": m, "split": s, **v} for m in mets for s, v in mets[m].items()]
table = pd.DataFrame(rows)
print("\n", table.round(3).to_string(index=False))
table.to_csv(f"{args.prefix}_metrics.csv", index=False)
pd.concat([
pd.DataFrame({"sample_idx": sets[s], "split": s, "observed": ys[s],
"pred_RF": preds["RF"][s], "pred_MLP": preds["MLP"][s]})
for s in sets
]).to_csv(f"{args.prefix}_predictions.csv", index=False)
make_figure(f"{args.prefix}_performance.png", preds, ys, mets, rf, wl, mlp)
print(f"\nSalvati: {args.prefix}_performance.png, {args.prefix}_metrics.csv, {args.prefix}_predictions.csv")
if __name__ == "__main__":
main()
Plottando i dati si vedono alcuni correlazioni (effettuate sul log) strette tra parametri
pH in acqua e CaCl2 sono praticamente identici e ridondanti
Carbonio Organico ed Azoto tot correlano a 0.98. La CEC dipende soprattutto da argilla e carbonio organico. Il CaCO3 è legato positivamente al pH, ma in modo non lineare: il pH satura sopra gli 8 circa. Il potassio estraibile tende a crescere con l’argilla. 
I dati sono stati elaborati tramite Random Tree Forest ed una tradizionale MLP (Multilayer Perceptron uno strato nascosto da 128 neuroni, due strati, da 256 e 64 neuroni e tre strati, da 512, 128 e 32 neuroni) Come si vede MLP funziona molto bene quando deve trasferire la conoscenza acquisita in train verso validation e test dataset (rispettivamente 70%,10% e 20% del dataset completo) model splitn R2 RMSE MAE RPIQ
RF train 28545 0.826 5.467 1.969 0.274
RF val 4078 0.625 8.129 3.310 0.185
RF test 8156 0.619 7.972 3.261 0.201
MLP train 28545 0.919 3.722 1.360 0.403
MLP val 4078 0.904 4.125 1.552 0.364
MLP test 8156 0.905 3.974 1.478 0.403
|
|---|
