Prosegue la serie degli schiaffi
 |
| Mappa random forest contro verita' a terra |
Ho provato a cambiare approccio passando random forest.
Per fare questa cosa pero' per avere un classificatore random forest oltre ad avere un dataset di positivi devo avere anche un dataset negativo
Per fare questo e' stato creato un tema puntuale di punti random non coincidenti con quelli del censimento amianto
import numpy as np
import geopandas as gpd
import shapely
from shapely.geometry import box
from scipy.spatial import cKDTree
# ---------------- parametri ----------------
SRC = "amianto_PRA_pubblico.shp" # punti amianto
MASK_RASTER = None # es. "mask_built.tif" (1 = edificato valido, stessa griglia EnMAP)
AREA = None # es. "area_campionamento.shp" (poligoni); ignorato se MASK_RASTER è impostato
MIN_DIST = 50 # metri minimi dai punti amianto
N_RATIO = 4 # negativi per ogni positivo
SEED = 0
OUT = "background_random.shp"
CRS_DEFAULT = 32632
# --------------------------------------------
rng = np.random.default_rng(SEED)
# punti amianto
pts = gpd.read_file(SRC)
pts = pts[pts.geometry.notna() & ~pts.geometry.is_empty]
if pts.crs is None:
pts = pts.set_crs(CRS_DEFAULT)
elif pts.crs.is_geographic:
pts = pts.to_crs(CRS_DEFAULT)
xy = np.column_stack([pts.geometry.x, pts.geometry.y])
n = len(xy)
n_bg = N_RATIO * n
tree = cKDTree(xy)
if MASK_RASTER:
# ---- campionamento su pixel della maschera (centri pixel, senza duplicati) ----
import rasterio
with rasterio.open(MASK_RASTER) as src:
if src.crs != pts.crs:
raise ValueError(f"CRS diverso: raster {src.crs} vs punti {pts.crs}. Riproietta uno dei due.")
mask = src.read(1) == 1
T = src.transform
rows, cols = np.nonzero(mask)
pick = rng.permutation(len(rows))
cx, cy = rasterio.transform.xy(T, rows[pick], cols[pick], offset="center")
cand = np.column_stack([cx, cy])
d, _ = tree.query(cand, k=1)
cand = cand[d > MIN_DIST]
if len(cand) < n_bg:
print(f"ATTENZIONE: disponibili solo {len(cand)} punti su {n_bg} richiesti. "
f"Riduci MIN_DIST o N_RATIO, oppure allarga la maschera.")
out = cand[:n_bg]
else:
# ---- campionamento casuale continuo dentro un poligono / bounding box ----
if AREA:
poly = gpd.read_file(AREA).to_crs(pts.crs).union_all()
else:
poly = box(*pts.total_bounds)
minx, miny, maxx, maxy = poly.bounds
out = np.empty((0, 2))
for _ in range(200): # limite di sicurezza contro cicli infiniti
if len(out) >= n_bg:
break
m = max(2 * (n_bg - len(out)), 10_000)
cand = np.column_stack([rng.uniform(minx, maxx, m),
rng.uniform(miny, maxy, m)])
cand = cand[shapely.contains_xy(poly, cand[:, 0], cand[:, 1])]
d, _ = tree.query(cand, k=1)
out = np.vstack([out, cand[d > MIN_DIST]])
if len(out) < n_bg:
print(f"ATTENZIONE: generati solo {len(out)} punti su {n_bg} richiesti. "
f"Riduci MIN_DIST o allarga l'area.")
out = out[:n_bg]
bg = gpd.GeoDataFrame(
{"id": np.arange(len(out)), "classe": 0},
geometry=gpd.points_from_xy(out[:, 0], out[:, 1]),
crs=pts.crs,
)
bg.to_file(OUT)
print(f"Positivi: {n} | Negativi salvati: {len(bg)} | file: {OUT}")
fatto questo sono stati estratti gli spettri dei rispettivi punti e creato un campo boolean false in caso di assenza di amianto. La libreria spettrale e' stata fusa con quella dei tetti in amianto ...quindi ho ottenuto una libreria con circa 14000 spettro con circa 50% presenza di amianto, 50% assenza di amianto
Questa libreria e' stata usata per addestrare una rete random forest con
import json, pickle
import numpy as np
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)
import joblib
# ---------------- parametri ----------------
GPKG = "C:/Users/l.innocenti/Documents/MF/random_forest.gpkg"
LAYER = "random_forest" # None = prima tabella; altrimenti nome tabella
FIELD_PROFILES = "profile"
FIELD_BG = "background"
CELL = 3000 # metri, cella per i gruppi della CV spaziale
# intervalli (nm) da escludere: vapore acqueo e bordi rumorosi
DROP_RANGES = [(0, 550), (1300, 1500), (1750, 2150), (2400, 3000)]
N_TREES = 500
SEED = 0
OUT_MODEL = "C:/Users/l.innocenti/Documents/MF/rf_amianto.joblib"
# --------------------------------------------
def parse_profile(v):
"""Decodifica il campo profile di EnMAP-Box/QPS (dict, JSON testo o blob)."""
if v is None:
return None
if isinstance(v, dict): # già decodificato (campo JSON del gpkg)
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
df = gpd.read_file(GPKG, layer=LAYER)
print("Record letti:", len(df), "| campi:", list(df.columns))
# --- spettri ---
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) # lunghezza più frequente
ok = np.array([l == n_bands for l in lengths])
print(f"Bande: {n_bands} | profili scartati per lunghezza/vuoti: {(~ok).sum()}")
first = next(p for p, o in zip(prof, ok) if o)
wl = np.array(first["x"], float) if first.get("x") is not None else None
bbl = np.array(first["bbl"], bool) if first.get("bbl") is not None else np.ones(n_bands, bool)
band_mask = bbl.copy()
if wl is not None:
for lo, hi in DROP_RANGES:
band_mask &= ~((wl >= lo) & (wl <= hi))
else:
print("ATTENZIONE: lunghezze d'onda assenti, uso solo bbl (nessun filtro per intervalli)")
X = np.array([np.asarray(p["y"], float) for p, o in zip(prof, ok) if o])
keep = np.where(ok)[0]
# --- etichette: 1 = amianto, 0 = background ---
bg = df[FIELD_BG].iloc[keep]
valid = bg.notna().to_numpy()
y = (~bg.fillna(False).astype(bool)).to_numpy().astype(int)
# --- pulizia spettri ---
X = X[:, band_mask]
good = valid & np.isfinite(X).all(axis=1) & (np.abs(X).sum(axis=1) > 0)
X, y, keep = X[good], y[good], keep[good]
print(f"Campioni: {len(y)} | amianto: {y.sum()} | background: {(y == 0).sum()} | bande usate: {X.shape[1]}")
# --- gruppi spaziali per la CV ---
groups = np.arange(len(y)) # fallback: ogni campione un gruppo
if hasattr(df, "geometry") and df.geometry.notna().any():
g = df.iloc[keep]
if g.crs is not None and g.crs.is_geographic:
g = g.to_crs(32632)
has = g.geometry.notna().to_numpy()
cx = (g.geometry.x.to_numpy()[has] // CELL).astype(np.int64)
cy = (g.geometry.y.to_numpy()[has] // CELL).astype(np.int64)
groups[has] = -(cx * 100000 + cy) - 1 # id negativi, distinti dal fallback
else:
print("ATTENZIONE: nessuna geometria, la CV NON è spaziale e le metriche saranno ottimistiche")
# --- modello ---
def make_rf():
return RandomForestClassifier(
n_estimators=N_TREES, max_features="sqrt", min_samples_leaf=3,
class_weight="balanced_subsample", n_jobs=-1, random_state=SEED)
# --- cross-validation spaziale ---
cv = StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=SEED)
oof = np.zeros(len(y))
for k, (tr, te) in enumerate(cv.split(X, y, groups), 1):
m = make_rf().fit(X[tr], y[tr])
oof[te] = m.predict_proba(X[te])[:, 1]
print(f"fold {k}: PR-AUC = {average_precision_score(y[te], oof[te]):.3f}")
print(f"\nPR-AUC globale: {average_precision_score(y, oof):.3f} | ROC-AUC: {roc_auc_score(y, oof):.3f}")
prec, rec, thr = precision_recall_curve(y, oof)
f1 = 2 * prec[:-1] * rec[:-1] / np.clip(prec[:-1] + rec[:-1], 1e-9, None)
i = f1.argmax()
print(f"Soglia F1 ottimale: {thr[i]:.3f} (precision {prec[i]:.3f}, recall {rec[i]:.3f})")
for t in (0.5, 0.7, 0.9):
j = np.searchsorted(thr, t)
print(f"soglia {t}: precision {prec[j]:.3f}, recall {rec[j]:.3f}")
# --- modello finale su tutti i dati ---
rf = make_rf().fit(X, y)
joblib.dump({"model": rf, "band_mask": band_mask, "wavelengths": wl,
"n_bands_total": n_bands, "threshold_f1": float(thr[i])}, OUT_MODEL)
print("Modello salvato in", OUT_MODEL)
# --- bande più importanti ---
imp = rf.feature_importances_
labels = wl[band_mask] if wl is not None else np.arange(X.shape[1])
for j in np.argsort(imp)[::-1][:10]:
print(f"banda {labels[j]:.0f}: {imp[j]:.4f}")
a questo punto ho preso una immagine sulla quale la rete non e' stata addestrata e fatto inferenza tramite
import numpy as np
import joblib
import rasterio
import geopandas as gpd
import xml.etree.ElementTree as ET
from rasterio.windows import Window
from scipy.ndimage import maximum_filter
from scipy.spatial import cKDTree
# ---------------- parametri ----------------
IMG = "C:/Users/l.innocenti/Documents/MF/ENMAP_SPECTRAL_IMAGE.TIF" # stesso prodotto da cui vengono gli spettri
METADATA_XML = "C:/Users/l.innocenti/Documents/MF/ENMAP_METADATA.XML" # METADATA.XML dello stesso prodotto
MODEL = "C:/Users/l.innocenti/Documents/MF/rf_amianto.joblib"
OUT_PROB = "C:/Users/l.innocenti/Documents/MF/prob_amianto.tif"
CENSUS = "amianto_PRA_pubblico.shp" # punti del censimento
LIBRARY = "C:/Users/l.innocenti/Documents/MF/random_forest.gpkg" # per escludere i punti usati nel training
LIBRARY_LAYER = "random_forest"
BUILT_MASK = None # raster 0/1 sulla stessa griglia dell'immagine (1 = edificato): rende il "caso" onesto
MTMF_SCORE = None # raster punteggio MTMF sulla stessa griglia, per il confronto a parità di pixel
NODATA_DEFAULT = -32768
BLOCK = 256
TOL_NM = 1.0 # scarto massimo ammesso tra lunghezza d'onda del modello e della banda abbinata
NEIGH = 2 # un punto è "hit" se c'è un pixel segnalato entro 2 pixel
EXCL_DIST = 150 # m: esclude i punti del censimento vicini a spettri usati nel training
TOP_FRACTIONS = [0.001, 0.005, 0.01, 0.02, 0.05] # quota di pixel validi segnalati
# --------------------------------------------
def wavelengths_from_xml(xml_path, n):
"""Legge un solo blocco da n lunghezze d'onda (nm) dal METADATA.XML di EnMAP."""
root = ET.parse(xml_path).getroot()
for tag in ("wavelengthCenterOfBand", "waveLength"):
arr = np.array([float(e.text) for e in root.iter()
if e.tag.split("}")[-1] == tag and e.text])
if len(arr) >= n and len(arr) % n == 0:
arr = arr[:n]
if 400 < arr[0] < 440 and 2400 < arr[-1] < 2500:
print(f"Lunghezze d'onda lette da <{tag}>: {arr[0]:.1f} ... {arr[-1]:.1f} nm")
return arr
raise SystemExit("Blocco di lunghezze d'onda non riconosciuto nell'XML")
# ===== 1) predizione a blocchi =====
d = joblib.load(MODEL)
rf, band_mask, n_model = d["model"], d["band_mask"], d["n_bands_total"]
wl_model = d["wavelengths"]
if wl_model is None:
raise SystemExit("Il modello non ha lunghezze d'onda salvate: impossibile abbinare le bande.")
with rasterio.open(IMG) as src:
wl_img = wavelengths_from_xml(METADATA_XML, src.count)
# abbinamento banda-modello -> banda-immagine, con indice crescente
# (il modello è un sottoinsieme ordinato delle bande dell'immagine)
img_idx = np.empty(len(wl_model), dtype=int)
prev = -1
for k, w in enumerate(wl_model):
cand = np.arange(prev + 1, len(wl_img))
if len(cand) == 0:
raise SystemExit("Abbinamento impossibile: finite le bande dell'immagine.")
j = cand[np.abs(wl_img[cand] - w).argmin()]
img_idx[k] = j
prev = j
err = np.abs(wl_img[img_idx] - wl_model)
print(f"Bande modello: {n_model} | bande immagine: {src.count} | "
f"scarto max: {err.max():.4f} nm | medio: {err.mean():.4f} nm")
if err.max() > TOL_NM:
raise SystemExit(f"Scarto oltre {TOL_NM} nm: le lunghezze d'onda non corrispondono, "
"libreria e immagine vengono da sensori/versioni diverse.")
if len(np.unique(img_idx)) != len(img_idx):
raise SystemExit("Due bande del modello cadono sulla stessa banda dell'immagine.")
missing = np.setdiff1d(np.arange(len(wl_img)), img_idx)
print(f"Bande dell'immagine non usate ({len(missing)}):", np.round(wl_img[missing], 1))
read_idx = (img_idx + 1).tolist() # rasterio: indici da 1
nodata = src.nodata if src.nodata is not None else NODATA_DEFAULT
H, W = src.height, src.width
T, CRS = src.transform, src.crs
prof = src.profile.copy()
prof.update(count=1, dtype="float32", nodata=-1.0, compress="lzw",
tiled=True, blockxsize=256, blockysize=256)
prob = np.full((H, W), -1.0, dtype=np.float32)
with rasterio.open(OUT_PROB, "w", **prof) as dst:
for r0 in range(0, H, BLOCK):
for c0 in range(0, W, BLOCK):
h, w = min(BLOCK, H - r0), min(BLOCK, W - c0)
win = Window(c0, r0, w, h)
arr = src.read(read_idx, window=win) # (n_model, h, w)
X = arr.reshape(len(read_idx), -1).T[:, band_mask].astype(np.float32)
good = (np.isfinite(X).all(axis=1)
& (X != nodata).all(axis=1)
& (np.abs(X).sum(axis=1) > 0))
p = np.full(X.shape[0], -1.0, dtype=np.float32)
if good.any():
p[good] = rf.predict_proba(X[good])[:, 1]
blk = p.reshape(h, w)
prob[r0:r0 + h, c0:c0 + w] = blk
dst.write(blk, 1, window=win)
print(f"righe {min(r0 + BLOCK, H)}/{H}")
print("Probabilità salvata in", OUT_PROB)
# ===== 2) confronto con il censimento =====
valid = prob >= 0
eval_mask = valid.copy()
if BUILT_MASK:
with rasterio.open(BUILT_MASK) as m:
built = m.read(1) == 1
if built.shape != valid.shape:
raise SystemExit("BUILT_MASK non ha la stessa griglia dell'immagine")
eval_mask &= built
else:
print("ATTENZIONE: senza BUILT_MASK il 'caso' è calcolato su tutti i pixel validi "
"(campagna, boschi...), quindi gli arricchimenti assoluti sono gonfiati")
cen = gpd.read_file(CENSUS)
cen = cen[cen.geometry.notna()].to_crs(CRS)
xs, ys = cen.geometry.x.to_numpy(), cen.geometry.y.to_numpy()
# esclude i punti vicini a spettri della libreria (positivi e background usati nel training)
lib = gpd.read_file(LIBRARY, layer=LIBRARY_LAYER)
lib = lib[lib.geometry.notna()].to_crs(CRS)
tree = cKDTree(np.column_stack([lib.geometry.x, lib.geometry.y]))
dist, _ = tree.query(np.column_stack([xs, ys]), k=1)
far = dist > EXCL_DIST
cols_f, rows_f = ~T * (xs, ys)
rows, cols = np.floor(rows_f).astype(int), np.floor(cols_f).astype(int)
inside = (rows >= 0) & (rows < H) & (cols >= 0) & (cols < W)
sel = inside & far
sel[sel] = eval_mask[rows[sel], cols[sel]] # solo pixel valutabili
rows, cols = rows[sel], cols[sel]
n_pts = len(rows)
print(f"\nPunti del censimento: {len(cen)} | nell'immagine: {inside.sum()} "
f"| lontani dal training e su pixel validi: {n_pts}")
def hit_stats(score, frac):
thr = np.quantile(score[eval_mask], 1 - frac)
flag = (score >= thr) & eval_mask
dil = maximum_filter(flag.astype(np.uint8), size=2 * NEIGH + 1) > 0
obs = dil[rows, cols].mean()
chance = dil[eval_mask].mean()
ci = 1.96 * np.sqrt(obs * (1 - obs) / n_pts)
return thr, int(flag.sum()), obs, ci, chance
scores = {"RF": prob}
if MTMF_SCORE:
with rasterio.open(MTMF_SCORE) as m:
s = m.read(1).astype(np.float32)
if s.shape != valid.shape:
raise SystemExit("MTMF_SCORE non ha la stessa griglia dell'immagine")
s[~np.isfinite(s)] = -np.inf
scores["MTMF"] = s
print(f"\n{'metodo':6} {'% pixel':>8} {'n pixel':>9} {'soglia':>9} {'hit %':>14} {'caso %':>8} {'arricch.':>9}")
for name, sc in scores.items():
for f in TOP_FRACTIONS:
thr, npx, obs, ci, chance = hit_stats(sc, f)
print(f"{name:6} {100*f:8.2f} {npx:9d} {thr:9.3f} "
f"{100*obs:6.1f}±{100*ci:4.1f} {100*chance:8.1f} {obs/max(chance,1e-9):9.2f}")
Punti del censimento: 56865 | nell'immagine: 4189 | lontani dal training e su pixel validi: 2925
metodo % pixel n pixel soglia hit % caso % arricch.
RF 0.10 1048 0.994 3.8± 0.7 1.2 3.23
RF 0.50 5238 0.970 14.9± 1.3 4.1 3.59
RF 1.00 10476 0.933 27.1± 1.6 6.7 4.03
RF 2.00 20952 0.876 47.4± 1.8 11.3 4.20
RF 5.00 52380 0.738 68.4± 1.7 21.3 3.22
Nella riga 2,00 segnali 20952 pixel, e il 47,4% dei punti del censimento ha un pixel segnalato a meno di 60 m. Se li avessi segnalati a caso, ti aspetteresti l'11,3%. Il rapporto è 47,4 / 11,3 ≈ 4,2.
Prendendo il file prob_amianto.tif e guardando quanti punti di verita' a terra cadono nell'intorno dei pixel
import numpy as np
import rasterio
import geopandas as gpd
from scipy.spatial import cKDTree
# ---------------- parametri ----------------
PROB = "C:/Users/l.innocenti/Documents/MF/prob_amianto.tif" # probabilità 0-1, nodata = -1
CENSUS = "C:/Users/l.innocenti/Documents/MF/amianto_PRA_pubblico.shp"
OUT = "C:/Users/l.innocenti/Documents/MF/punti_distanza.gpkg" # None per non salvare
DIST = 60.0 # metri
THRESHOLDS = [0.7, 0.8, 0.9, 0.95, 0.99] # un pixel è "segnalato" se prob >= soglia
N_CHANCE = 200_000 # pixel validi campionati per stimare il "caso"
SEED = 0
# opzionale: escludi i punti entro EXCL_DIST m da spettri usati nel training
LIBRARY = None # es. "C:/.../random_forest.gpkg"
LIBRARY_LAYER = "random_forest"
EXCL_DIST = 150.0
# --------------------------------------------
with rasterio.open(PROB) as src:
prob = src.read(1)
T, CRS = src.transform, src.crs
H, W = src.height, src.width
px = abs(T.a)
nodata = src.nodata if src.nodata is not None else -1.0
valid = np.isfinite(prob) & (prob != nodata) & (prob >= 0)
print(f"Immagine {W}x{H} | pixel {px:.0f} m | pixel validi: {valid.sum()}")
# ---------- punti ----------
cen = gpd.read_file(CENSUS)
cen = cen[cen.geometry.notna()].to_crs(CRS).reset_index(drop=True)
xs, ys = cen.geometry.x.to_numpy(), cen.geometry.y.to_numpy()
n_tot = len(cen)
cf, rf_ = ~T * (xs, ys)
cols, rows = np.floor(cf).astype(int), np.floor(rf_).astype(int)
inside = (rows >= 0) & (rows < H) & (cols >= 0) & (cols < W)
use = inside.copy()
use[inside] = valid[rows[inside], cols[inside]] # solo punti su pixel validi
print(f"Punti totali: {n_tot} | dentro l'immagine: {inside.sum()} | su pixel validi: {use.sum()}")
if LIBRARY:
lib = gpd.read_file(LIBRARY, layer=LIBRARY_LAYER)
lib = lib[lib.geometry.notna()].to_crs(CRS)
d_lib, _ = cKDTree(np.column_stack([lib.geometry.x, lib.geometry.y])).query(
np.column_stack([xs, ys]), k=1)
use &= d_lib > EXCL_DIST
print(f"Dopo esclusione punti vicini al training (< {EXCL_DIST:.0f} m): {use.sum()}")
pts = np.column_stack([xs[use], ys[use]])
n = len(pts)
if n == 0:
raise SystemExit("Nessun punto utilizzabile.")
# ---------- "caso": pixel validi casuali ----------
rng = np.random.default_rng(SEED)
vr, vc = np.nonzero(valid)
pick = rng.choice(len(vr), size=min(N_CHANCE, len(vr)), replace=False)
cx, cy = T * (vc[pick] + 0.5, vr[pick] + 0.5)
rand_pts = np.column_stack([cx, cy])
# ---------- distanza dal pixel segnalato più vicino (centro pixel) ----------
res = cen.loc[use, ["geometry"]].copy()
print(f"\nDistanza dal centro del pixel segnalato più vicino, soglia {DIST:.0f} m "
f"(n punti = {n})\n")
print(f"{'soglia':>7} {'pixel segn.':>12} {'< %dm' % DIST:>8} {'>= %dm' % DIST:>8} "
f"{'% entro':>8} {'caso %':>8} {'arricch.':>9}")
for t in THRESHOLDS:
fr, fc = np.nonzero(valid & (prob >= t))
if len(fr) == 0:
print(f"{t:7.2f} {0:12d} (nessun pixel sopra soglia)")
continue
fx, fy = T * (fc + 0.5, fr + 0.5)
tree = cKDTree(np.column_stack([fx, fy]))
d, _ = tree.query(pts, k=1)
d_rand, _ = tree.query(rand_pts, k=1)
near = d < DIST
obs = near.mean()
chance = (d_rand < DIST).mean()
print(f"{t:7.2f} {len(fr):12d} {near.sum():8d} {(~near).sum():8d} "
f"{100 * obs:8.1f} {100 * chance:8.1f} {obs / max(chance, 1e-9):9.2f}")
res[f"d_{t:.2f}"] = d
# ---------- salvataggio ----------
if OUT:
res.to_file(OUT, driver="GPKG")
print("\nSalvato:", OUT, "(distanza in metri per ogni soglia)")
si ha che con una soglia di 0.7 si ha una percentuale di corretta detection entro 60 m del 65% dei casi (tirando a caso la percentuale sarebbe stata del 14.5% quindi il segnale e' ben presente)
Rendendo la soglia piu' stringente la percentuale crolla miseramente
soglia pixel segn. < 60m >= 60m % entro caso % arricch.
0.70 60332 1912 1013 65.4 14.5 4.51
0.80 38932 1590 1335 54.4 10.3 5.29
0.90 16122 959 1966 32.8 5.2 6.28
0.95 8183 514 2411 17.6 3.0 5.82
0.99 1747 109 2816 3.7 0.8 4.43
