Per cercare di migliorare le librerie dell'amianto su Enmap ho provato ad essere piu' selettivo
Usando WorldCover Esa Map ho selezionato in QGis i build-up che sono indicati da Banda 1 Tavolozza 50
ho tagliato il raster usando come base l'impronta dell'immagine Enmap con Processing → Toolbox → GDAL → Raster extraction → Clip raster by extent
alla fine ho ottenuto un raster 0-1 in cui 1 indica la posizione dei built-up e 0 tutte le altre classiche
A questo punto ho convertito il raster in un tema puntuale tramite Processing → Toolbox → GDAL → Conversione raster → Pixel raster in punti
Alla fine ho un tema puntuale di tetti..a questo punto devo selezionare solo i tetti che non sono di amianto
Per questo motivo ho creato un buffer di 60 metri attorno ai punti del tema con i punti dei tetti di amianto_PRA_pubblico con Processing → Toolbox → Geometria vettoriale → Buffer
siamo quasi alla fine, togliendo dal tema tutti i pixel che sono esterni al buffer si puo' usare Processing → Toolbox → Estrai per posizione
Per estrarre dei punti casuali all'interno del tema ho usato Processing Toolobox Estrazione del vettore Estrazione casuale estraendo un nuovo tema di 20000 punti
In seguito Enmap Toolbox Spectral Library Extract spectral profiles from raster layer per avere una libreria spettrale di tetti non in amianto
Come nel precedente post ho creato un campo boolean amianto (true per i tetti in amianto, false l'ovvio opposto) ed ho unito le due librerie con
Per unire la libreria dei tetti in amianto e quelle non in amianto si apre la libreria A, si selezionano tutti i profili si copiano in Spectral Library Viewer di Enmap Toolbox e si incollano nella libreria B
Per rendere il modello misurabile la libreria e' stata splittata in 70% train, 20% test e 10% validate
import json, pickle
import numpy as np
import pandas as pd
import geopandas as gpd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.metrics import (average_precision_score, roc_auc_score,
precision_recall_curve, precision_score,
recall_score, f1_score)
# ---------------- parametri ----------------
GPKG = "C:/Users/l.innocenti/Documents/MF/worldview/libreria/unione.gpkg"
LAYER = "no_amianto2"
FIELD_PROFILES = "profile"
FIELD_BG = "amianto" # 1.0 = amianto; NaN/vuoto = background
CELL = 500 # metri
EDGE_DROP = [(0, 550), (2400, 3000)] # sempre escluse
WATER_DROP = [(1300, 1500), (1750, 2150)] # escluse solo nelle varianti "senza bande"
N_TREES = 500
SEEDS = [0, 1, 2] # ogni seed cambia sia lo split spaziale sia il random forest
N_FOLDS = 10
FOLDS_TRAIN, FOLDS_TEST, FOLDS_VAL = (0, 1, 2, 3, 4, 5, 6), (7, 8), (9,)
OUT_CSV = "C:/Users/l.innocenti/Documents/MF/ablazione_risultati.csv"
TRUE_VALUES = {"true", "1", "1.0", "yes", "si", "sì", "t"}
# --------------------------------------------
def parse_profile(v):
if v is None:
return None
if isinstance(v, dict):
return v
if isinstance(v, (bytes, bytearray, memoryview)):
b = bytes(v)
try:
return json.loads(b.decode("utf-8"))
except Exception:
return pickle.loads(b) # solo su file tuoi, di cui ti fidi
if isinstance(v, str):
return json.loads(v)
return None
def normalize_spectra(A, mode):
if mode is None:
return A
n = np.linalg.norm(A, axis=1, keepdims=True)
return A / np.clip(n, 1e-12, None)
def make_rf(seed):
return RandomForestClassifier(
n_estimators=N_TREES, max_features="sqrt", min_samples_leaf=3,
class_weight="balanced_subsample", n_jobs=-1, random_state=seed)
# --- lettura ---
df = gpd.read_file(GPKG, layer=LAYER)
print("Record letti:", len(df), "| campi:", list(df.columns))
print(df[FIELD_BG].value_counts(dropna=False))
prof = [parse_profile(v) for v in df[FIELD_PROFILES]]
lengths = [len(p["y"]) if p and p.get("y") is not None else 0 for p in prof]
n_bands = max(set(lengths), key=lengths.count)
ok = np.array([l == n_bands for l in lengths])
first = next(p for p, o in zip(prof, ok) if o)
wl = np.array(first["x"], float)
bbl = np.array(first["bbl"], bool) if first.get("bbl") is not None else np.ones(n_bands, bool)
# maschera larga (con bande reintrodotte) e stretta (senza)
mask_wide = bbl.copy()
for lo, hi in EDGE_DROP:
mask_wide &= ~((wl >= lo) & (wl <= hi))
mask_narrow = mask_wide.copy()
for lo, hi in WATER_DROP:
mask_narrow &= ~((wl >= lo) & (wl <= hi))
print(f"Bande: larga {mask_wide.sum()} | stretta {mask_narrow.sum()}")
X = np.array([np.asarray(p["y"], float) for p, o in zip(prof, ok) if o])
keep = np.where(ok)[0]
lab = df[FIELD_BG].iloc[keep]
y = lab.map(lambda v: 1 if str(v).strip().lower() in TRUE_VALUES else 0).to_numpy().astype(int)
# stessi campioni per tutte le varianti: filtro sulla maschera larga
Xw = X[:, mask_wide]
good = np.isfinite(Xw).all(axis=1) & (np.abs(Xw).sum(axis=1) > 0)
X, y, keep, Xw = X[good], y[good], keep[good], Xw[good]
assert len(Xw) == len(y) == len(keep), "Spettri e etichette non allineati"
print(f"Campioni: {len(y)} | amianto: {y.sum()} | background: {(y == 0).sum()}")
if y.sum() == 0 or (y == 0).sum() == 0:
raise SystemExit("Una sola classe presente: controlla FIELD_BG")
# indici delle colonne della maschera stretta dentro la larga
sub_narrow = mask_narrow[mask_wide]
# --- gruppi spaziali ---
g = df.iloc[keep]
if g.crs is not None and g.crs.is_geographic:
g = g.to_crs(32632)
cx = (g.geometry.x.to_numpy() // CELL).astype(np.int64)
cy = (g.geometry.y.to_numpy() // CELL).astype(np.int64)
groups = -(cx * 100000 + cy) - 1
print("Celle spaziali:", len(np.unique(groups)))
configs = [("grezzo", None, False), ("grezzo", None, True),
("L2", "l2", False), ("L2", "l2", True)]
rows = []
for seed in SEEDS:
sgkf = StratifiedGroupKFold(n_splits=N_FOLDS, shuffle=True, random_state=seed)
fold_id = np.zeros(len(y), int)
for k, (_, te) in enumerate(sgkf.split(X, y, groups)):
fold_id[te] = k
tr = np.where(np.isin(fold_id, FOLDS_TRAIN))[0]
te = np.where(np.isin(fold_id, FOLDS_TEST))[0]
va = np.where(np.isin(fold_id, FOLDS_VAL))[0]
for norm_name, norm, wide in configs:
Xc = Xw if wide else Xw[:, sub_narrow]
Xc = normalize_spectra(Xc, norm)
rf = make_rf(seed).fit(Xc[tr], y[tr])
pv = rf.predict_proba(Xc[va])[:, 1]
pt = rf.predict_proba(Xc[te])[:, 1]
prec, rec, thr = precision_recall_curve(y[va], pv)
f1c = 2 * prec[:-1] * rec[:-1] / np.clip(prec[:-1] + rec[:-1], 1e-9, None)
t = float(thr[f1c.argmax()])
yp = (pt >= t).astype(int)
rows.append({
"seed": seed, "spettri": norm_name,
"bande_1300-1500/1750-2150": "incluse" if wide else "escluse",
"n_bande": Xc.shape[1],
"PR-AUC": average_precision_score(y[te], pt),
"ROC-AUC": roc_auc_score(y[te], pt),
"soglia": t,
"precision": precision_score(y[te], yp, zero_division=0),
"recall": recall_score(y[te], yp),
"F1": f1_score(y[te], yp),
})
r = rows[-1]
print(f"seed {seed} | {norm_name:6s} | bande {r['bande_1300-1500/1750-2150']:7s} "
f"| PR-AUC {r['PR-AUC']:.3f} | ROC-AUC {r['ROC-AUC']:.3f} | F1 {r['F1']:.3f}")
res = pd.DataFrame(rows)
res.to_csv(OUT_CSV, index=False)
keys = ["spettri", "bande_1300-1500/1750-2150", "n_bande"]
metrics = ["PR-AUC", "ROC-AUC", "precision", "recall", "F1"]
summary = res.groupby(keys)[metrics].agg(["mean", "std"]).round(3)
pd.set_option("display.width", 250)
pd.set_option("display.max_columns", 50)
print(f"\n=== RISULTATI SUL TEST (media e dev. std su {len(SEEDS)} seed, CELL={CELL} m) ===")
print(summary)
print(f"\nBaseline PR-AUC (prevalenza amianto): {y.mean():.3f}")
print("Risultati per singolo seed salvati in", OUT_CSV)
import json, pickle
import numpy as np
import pandas as pd
import geopandas as gpd
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.metrics import (average_precision_score, roc_auc_score,
precision_recall_curve, precision_score,
recall_score, f1_score)
# ---------------- parametri ----------------
GPKG = "C:/Users/l.innocenti/Documents/MF/worldview/libreria/unione.gpkg"
LAYER = "no_amianto2"
FIELD_PROFILES = "profile"
FIELD_BG = "amianto" # 1.0 = amianto; NaN/vuoto = background
CELL = 500 # metri
EDGE_DROP = [(0, 550), (2400, 3000)] # sempre escluse
WATER_DROP = [(1300, 1500), (1750, 2150)] # escluse solo nelle varianti "senza bande"
N_TREES = 500
SEEDS = [0, 1, 2] # ogni seed cambia sia lo split spaziale sia il random forest
N_FOLDS = 10
FOLDS_TRAIN, FOLDS_TEST, FOLDS_VAL = (0, 1, 2, 3, 4, 5, 6), (7, 8), (9,)
OUT_CSV = "C:/Users/l.innocenti/Documents/MF/ablazione_risultati.csv"
TRUE_VALUES = {"true", "1", "1.0", "yes", "si", "sì", "t"}
# --------------------------------------------
def parse_profile(v):
if v is None:
return None
if isinstance(v, dict):
return v
if isinstance(v, (bytes, bytearray, memoryview)):
b = bytes(v)
try:
return json.loads(b.decode("utf-8"))
except Exception:
return pickle.loads(b) # solo su file tuoi, di cui ti fidi
if isinstance(v, str):
return json.loads(v)
return None
def normalize_spectra(A, mode):
if mode is None:
return A
n = np.linalg.norm(A, axis=1, keepdims=True)
return A / np.clip(n, 1e-12, None)
def make_rf(seed):
return RandomForestClassifier(
n_estimators=N_TREES, max_features="sqrt", min_samples_leaf=3,
class_weight="balanced_subsample", n_jobs=-1, random_state=seed)
# --- lettura ---
df = gpd.read_file(GPKG, layer=LAYER)
print("Record letti:", len(df), "| campi:", list(df.columns))
print(df[FIELD_BG].value_counts(dropna=False))
prof = [parse_profile(v) for v in df[FIELD_PROFILES]]
lengths = [len(p["y"]) if p and p.get("y") is not None else 0 for p in prof]
n_bands = max(set(lengths), key=lengths.count)
ok = np.array([l == n_bands for l in lengths])
first = next(p for p, o in zip(prof, ok) if o)
wl = np.array(first["x"], float)
bbl = np.array(first["bbl"], bool) if first.get("bbl") is not None else np.ones(n_bands, bool)
# maschera larga (con bande reintrodotte) e stretta (senza)
mask_wide = bbl.copy()
for lo, hi in EDGE_DROP:
mask_wide &= ~((wl >= lo) & (wl <= hi))
mask_narrow = mask_wide.copy()
for lo, hi in WATER_DROP:
mask_narrow &= ~((wl >= lo) & (wl <= hi))
print(f"Bande: larga {mask_wide.sum()} | stretta {mask_narrow.sum()}")
X = np.array([np.asarray(p["y"], float) for p, o in zip(prof, ok) if o])
keep = np.where(ok)[0]
lab = df[FIELD_BG].iloc[keep]
y = lab.map(lambda v: 1 if str(v).strip().lower() in TRUE_VALUES else 0).to_numpy().astype(int)
# stessi campioni per tutte le varianti: filtro sulla maschera larga
Xw = X[:, mask_wide]
good = np.isfinite(Xw).all(axis=1) & (np.abs(Xw).sum(axis=1) > 0)
X, y, keep, Xw = X[good], y[good], keep[good], Xw[good]
assert len(Xw) == len(y) == len(keep), "Spettri e etichette non allineati"
print(f"Campioni: {len(y)} | amianto: {y.sum()} | background: {(y == 0).sum()}")
if y.sum() == 0 or (y == 0).sum() == 0:
raise SystemExit("Una sola classe presente: controlla FIELD_BG")
# indici delle colonne della maschera stretta dentro la larga
sub_narrow = mask_narrow[mask_wide]
# --- gruppi spaziali ---
g = df.iloc[keep]
if g.crs is not None and g.crs.is_geographic:
g = g.to_crs(32632)
cx = (g.geometry.x.to_numpy() // CELL).astype(np.int64)
cy = (g.geometry.y.to_numpy() // CELL).astype(np.int64)
groups = -(cx * 100000 + cy) - 1
print("Celle spaziali:", len(np.unique(groups)))
configs = [("grezzo", None, False), ("grezzo", None, True),
("L2", "l2", False), ("L2", "l2", True)]
rows = []
for seed in SEEDS:
sgkf = StratifiedGroupKFold(n_splits=N_FOLDS, shuffle=True, random_state=seed)
fold_id = np.zeros(len(y), int)
for k, (_, te) in enumerate(sgkf.split(X, y, groups)):
fold_id[te] = k
tr = np.where(np.isin(fold_id, FOLDS_TRAIN))[0]
te = np.where(np.isin(fold_id, FOLDS_TEST))[0]
va = np.where(np.isin(fold_id, FOLDS_VAL))[0]
for norm_name, norm, wide in configs:
Xc = Xw if wide else Xw[:, sub_narrow]
Xc = normalize_spectra(Xc, norm)
rf = make_rf(seed).fit(Xc[tr], y[tr])
pv = rf.predict_proba(Xc[va])[:, 1]
pt = rf.predict_proba(Xc[te])[:, 1]
prec, rec, thr = precision_recall_curve(y[va], pv)
f1c = 2 * prec[:-1] * rec[:-1] / np.clip(prec[:-1] + rec[:-1], 1e-9, None)
t = float(thr[f1c.argmax()])
yp = (pt >= t).astype(int)
rows.append({
"seed": seed, "spettri": norm_name,
"bande_1300-1500/1750-2150": "incluse" if wide else "escluse",
"n_bande": Xc.shape[1],
"PR-AUC": average_precision_score(y[te], pt),
"ROC-AUC": roc_auc_score(y[te], pt),
"soglia": t,
"precision": precision_score(y[te], yp, zero_division=0),
"recall": recall_score(y[te], yp),
"F1": f1_score(y[te], yp),
})
r = rows[-1]
print(f"seed {seed} | {norm_name:6s} | bande {r['bande_1300-1500/1750-2150']:7s} "
f"| PR-AUC {r['PR-AUC']:.3f} | ROC-AUC {r['ROC-AUC']:.3f} | F1 {r['F1']:.3f}")
res = pd.DataFrame(rows)
res.to_csv(OUT_CSV, index=False)
keys = ["spettri", "bande_1300-1500/1750-2150", "n_bande"]
metrics = ["PR-AUC", "ROC-AUC", "precision", "recall", "F1"]
summary = res.groupby(keys)[metrics].agg(["mean", "std"]).round(3)
pd.set_option("display.width", 250)
pd.set_option("display.max_columns", 50)
print(f"\n=== RISULTATI SUL TEST (media e dev. std su {len(SEEDS)} seed, CELL={CELL} m) ===")
print(summary)
print(f"\nBaseline PR-AUC (prevalenza amianto): {y.mean():.3f}")
print("Risultati per singolo seed salvati in", OUT_CSV)
In conclusione il modello ha trovato 642 tetti in amianto su 1014 (ne ha mancati 372) e ha segnalato per errore 806 pixel di background su 2542.
Precision = TP / (TP + FP): 0.73 (tra i pixel che il modello segnala come amianto, quanti lo sono davvero)
Recall = TP / (TP + FN): 0.057 (pixel che sono amianto, quanti il modello ne trova)
Bande più importanti:
banda 1449: 0.0112
banda 1128: 0.0099
banda 2225: 0.0089
banda 1140: 0.0085
banda 1996: 0.0083
banda 2216: 0.0081
banda 1152: 0.0079
banda 951: 0.0079
banda 941: 0.0079
banda 1259: 0.0078
Nessun commento:
Posta un commento