Ho trovato un database pubblico di tetti in amianto pubblicato dalla Regione Piemonte ed ho provato ad incrociarlo con i dati di Enmap
Il primo test e' stato quello di estrarre le firme spettrali corrispondenti ad ogni pixel (dopo aver filtrato con Savitzky Gorlay) usando Enmap Toolbox con Extract spectral profiles from raster
L'immagine e' stata convertita in formato Envi Bil con Save raster layer as con il profile Envi BIL
Visto che molti pixel non sono puri ho provato l'algoritmo MTMF che lavora sub pixel
Prendendo lo spettro medio della libreria ho provato a calcolare l'algoritmo su tutta l'immagine tramite lo script sottostante che produce una mappa di MF Score, una di Infeasibility, una mappa di sintesi delle due precedenti con score 0-100 ed una mappa a 3 classi
import numpy as np
import rasterio
from spectral.io import envi
NODATA_OUT = -9999.0
def load_target(lib_hdr_path):
lib = envi.open(lib_hdr_path)
sp = np.asarray(lib.spectra, dtype=np.float64)
ok = np.isfinite(sp).all(axis=1) & (sp > -1000).all(axis=1)
print(f"Spettri validi in libreria: {ok.sum()} su {sp.shape[0]}")
return sp[ok]
def estimate_noise_cov(cube, valid2d):
"""Covarianza del rumore con shift-difference orizzontale. cube: (h, w, nb)."""
d = cube[:, 1:, :] - cube[:, :-1, :]
m = valid2d[:, 1:] & valid2d[:, :-1]
d = d[m].astype(np.float64)
return np.cov(d, rowvar=False) / 2.0
def run_mtmf(
image_path,
lib_hdr_path,
out_path="mtmf_output_test2.bsq",
lib_scale=None, # None = auto (allinea la libreria all'immagine)
mnf_min_eig=2.0, # tieni componenti MNF con autovalore > soglia (SNR > ~1)
max_comp=40,
z_lo=3.0, z_hi=8.0, # MF: z-score robusto dove il termine MF va da 0 a 1
inf_p_lo=5, inf_p_hi=75, # percentili di infeasibility: <=p_lo -> 1, >=p_hi -> 0
class_thr=(25, 50, 75), # soglie score per classi 1/2/3
plot=True,
):
# ---------- Immagine ----------
with rasterio.open(image_path) as src:
img = src.read().astype(np.float32)
nb, h, w = img.shape
meta = src.meta.copy()
nodata = src.nodata
pixels = np.transpose(img, (1, 2, 0)).reshape(-1, nb)
del img
valid = np.isfinite(pixels).all(axis=1)
if nodata is not None:
valid &= (pixels != nodata).all(axis=1)
valid &= (pixels != 0).any(axis=1)
print(f"Pixel validi: {valid.sum()} su {pixels.shape[0]}")
# ---------- Target ----------
spectra = load_target(lib_hdr_path)
if spectra.shape[1] != nb:
print("Numero di bande diverso tra immagine e libreria")
return
target = spectra.mean(axis=0)
vmax = np.nanmax(pixels[valid][:, ::10])
if lib_scale is None:
lib_scale = 10000.0 if (vmax > 2 and np.nanmax(target) <= 2) else 1.0
target = target * lib_scale
print(f"Fattore scala libreria: {lib_scale}")
# ---------- Bande utilizzabili ----------
vp = pixels[valid]
keep = (vp.std(axis=0) > 0) & np.isfinite(target)
nk = int(keep.sum())
print(f"Bande usate: {nk} su {nb}")
vp = vp[:, keep].astype(np.float64)
t = target[keep]
# ---------- Rumore + MNF ----------
cube_k = pixels[:, keep].reshape(h, w, nk)
Cn = estimate_noise_cov(cube_k, valid.reshape(h, w))
del cube_k, pixels
s, U = np.linalg.eigh(Cn)
s = np.clip(s, s.max() * 1e-8, None)
P = U / np.sqrt(s) # sbianca il rumore
mu = vp.mean(axis=0)
Xc = vp - mu
del vp
Cd = (Xc.T @ Xc) / (Xc.shape[0] - 1)
lam, V = np.linalg.eigh(P.T @ Cd @ P)
order = np.argsort(lam)[::-1]
lam, V = lam[order], V[:, order]
T = P @ V # trasformazione MNF completa
k = int(min(max_comp, (lam > mnf_min_eig).sum()))
k = max(k, 3)
print(f"Componenti MNF usate: {k} (autovalori: {np.round(lam[:k], 1)[:8]} ...)")
Y = Xc @ T[:, :k]
del Xc
d = (t - mu) @ T[:, :k] # target nello spazio MNF (centrato)
# ---------- Matched Filter ----------
Sigma = np.cov(Y, rowvar=False)
Sinv = np.linalg.pinv(Sigma)
a = Sinv @ d
denom = d @ a
if not np.isfinite(denom) or denom <= 0:
print("Denominatore non valido")
return
mf = (Y @ a) / denom # 0 = background, 1 = target
# ---------- Infeasibility (unita' di sigma di rumore) ----------
resid = Y - np.outer(mf, d)
infeas = np.sqrt((resid ** 2).sum(axis=1) / k)
del resid, Y
# ---------- Score di sintesi ----------
med = np.median(mf)
mad = 1.4826 * np.median(np.abs(mf - med))
z = (mf - med) / mad
mf_term = np.clip((z - z_lo) / (z_hi - z_lo), 0, 1)
lo, hi = np.percentile(infeas, [inf_p_lo, inf_p_hi])
feas_term = 1.0 - np.clip((infeas - lo) / (hi - lo), 0, 1)
score = 100.0 * mf_term * feas_term
classes = np.digitize(score, class_thr).astype(np.float32) # 0..3
print(f"MF: mediana={med:.3f}, sigma robusto={mad:.3f}")
print(f"Infeasibility: p{inf_p_lo}={lo:.2f}, p{inf_p_hi}={hi:.2f}")
for c, name in enumerate(["nessuna", "bassa", "media", "alta"]):
print(f" classe {c} ({name}): {(classes == c).sum()} pixel")
# ---------- Scrittura ----------
def to_full(arr):
out = np.full(h * w, NODATA_OUT, dtype=np.float32)
out[valid] = arr.astype(np.float32)
out[~np.isfinite(out)] = NODATA_OUT
return out.reshape(h, w)
meta.update(count=4, dtype="float32", driver="ENVI",
interleave="bsq", nodata=NODATA_OUT)
with rasterio.open(out_path, "w", **meta) as dst:
dst.write(to_full(mf), 1)
dst.write(to_full(infeas), 2)
dst.write(to_full(score), 3)
dst.write(to_full(classes), 4)
for i, n in enumerate(["MF score", "Infeasibility (sigma)",
"MTMF score 0-100", "MTMF classe 0-3"], 1):
dst.set_band_description(i, n)
print(f"Salvato: {out_path}")
# ---------- Scatter MF vs infeasibility ----------
if plot:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
idx = np.random.default_rng(0).choice(mf.size, min(200000, mf.size), replace=False)
plt.figure(figsize=(7, 5))
plt.scatter(mf[idx], infeas[idx], s=1, c=score[idx], cmap="viridis", alpha=0.5)
plt.colorbar(label="MTMF score")
plt.xlabel("MF score")
plt.ylabel("Infeasibility (sigma)")
plt.xlim(np.percentile(mf, 0.1), np.percentile(mf, 99.99))
plt.tight_layout()
plt.savefig("mtmf_scatter.png", dpi=150)
print("Salvato: mtmf_scatter.png")
except Exception as e:
print("Scatter non generato:", e)
if __name__ == "__main__":
run_mtmf("test2.bil", "libreria.hdr")
se si incrocia la mappa dello score sintetico con i punti verita' a terra tramite lo script sottostante ho una corrispondenza di circa il 60% di positivi correttamente individuati a terra
import numpy as np
import rasterio
import geopandas as gpd
from scipy.ndimage import distance_transform_edt
RASTER = "mtmf_output_test2.bsq"
SHAPE = "amianto_PRA_pubblico.shp"
SCORE_MIN = 6.0 # soglia: score > 6
MAX_DIST = 1.5 # distanza massima in pixel (strettamente < MAX_DIST)
BAND = 3 # MTMF score
with rasterio.open(RASTER) as src:
score = src.read(BAND).astype(np.float32)
nodata = src.nodata
crs = src.crs
h, w = score.shape
valid = np.isfinite(score)
if nodata is not None:
valid &= score != nodata
mask = valid & (score > SCORE_MIN)
print(f"Pixel con score > {SCORE_MIN}: {mask.sum()} su {valid.sum()} validi "
f"({100 * mask.sum() / valid.sum():.2f}%)")
if mask.sum() == 0:
raise SystemExit("Nessun pixel sopra soglia")
# distanza (in pixel) dal pixel "positivo" piu' vicino
dist = distance_transform_edt(~mask)
# ---------- Punti ----------
gdf = gpd.read_file(SHAPE)
print(f"Punti nello shapefile: {len(gdf)} (CRS: {gdf.crs})")
if gdf.crs != crs:
gdf = gdf.to_crs(crs)
with rasterio.open(RASTER) as src:
xs = gdf.geometry.x.values
ys = gdf.geometry.y.values
rows, cols = rasterio.transform.rowcol(src.transform, xs, ys)
rows = np.asarray(rows)
cols = np.asarray(cols)
inside = (rows >= 0) & (rows < h) & (cols >= 0) & (cols < w)
d = np.full(len(gdf), np.nan)
on_valid = np.zeros(len(gdf), dtype=bool)
d[inside] = dist[rows[inside], cols[inside]]
on_valid[inside] = valid[rows[inside], cols[inside]]
hit = inside & on_valid & (d < MAX_DIST)
gdf["dist_px"] = d
gdf["hit"] = hit.astype(int)
n_in = int(inside.sum())
n_val = int((inside & on_valid).sum())
print(f"\nPunti dentro l'immagine: {n_in}")
print(f"Punti su pixel validi (non nodata): {n_val}")
print(f"Punti a < {MAX_DIST:g} pixel da un pixel con score > {SCORE_MIN}: "
f"{hit.sum()} ({100 * hit.sum() / max(n_val, 1):.1f}% dei punti validi)")
# ---------- Baseline: quanto ci si aspetterebbe per caso ----------
area_frac = ((dist < MAX_DIST) & valid).sum() / valid.sum()
print(f"\nBaseline: {100 * area_frac:.1f}% dei pixel validi dell'immagine sta a "
f"< {MAX_DIST:g} pixel da un pixel sopra soglia")
print(f"Attesi per puro caso: ~{area_frac * n_val:.0f} punti su {n_val}")
# ---------- Salvataggio ----------
gdf.to_file("amianto_PRA_mtmf_check.gpkg", driver="GPKG")
print("Salvato: amianto_PRA_mtmf_check.gpkg (campi dist_px, hit)")
facendo una verifica con una immagine sulla quale la libreria non e' stata addestrata (la track vicina) la percentuale crolla ad misero 3% .
Dopo un po' di ricerca il problema e' chiaro: MTMF usa la covarianza interna all'immagine, se si cambia l'immagine cambia il background e quindi l'addestramento non e' piu' utile
Ho provato un approccio indipendente dall'immagine ovvero considerando la profondita' di picco tra 2320 e 2330 nm
#!/usr/bin/env python3
"""
Estrae il valore della profondita' di banda (raster prodotto da band_depth_2323.py)
nei punti di uno shapefile e calcola media e deviazione standard.
Uso
---
python extract_depth_at_points.py depth_2323.tif amianto_PRA_pubblico.shp
python extract_depth_at_points.py depth_wide.tif amianto_PRA_pubblico.shp --radius 1
python extract_depth_at_points.py depth_2323.tif punti.shp -o punti_depth.gpkg
NB: l'input raster e' il GeoTIFF di OUTPUT di band_depth_2323.py (es. depth_2323.tif),
non lo script .py.
Requisiti: numpy, rasterio, geopandas
"""
import argparse
import sys
import numpy as np
import rasterio
import geopandas as gpd
def neighborhood(arr, r, stat):
"""Statistica in una finestra quadrata (2r+1) ignorando i NaN.
Il valore e' assegnato solo dove il pixel centrale e' valido."""
if r == 0:
return arr
H, W = arr.shape
pad = np.pad(arr, r, constant_values=np.nan)
if stat == "max":
out = np.full(arr.shape, np.nan, dtype=np.float32)
for dy in range(2 * r + 1):
for dx in range(2 * r + 1):
out = np.fmax(out, pad[dy:dy + H, dx:dx + W])
else: # mean
s = np.zeros(arr.shape, dtype=np.float64)
n = np.zeros(arr.shape, dtype=np.float64)
for dy in range(2 * r + 1):
for dx in range(2 * r + 1):
w = pad[dy:dy + H, dx:dx + W]
ok = np.isfinite(w)
s[ok] += w[ok]
n[ok] += 1
with np.errstate(invalid="ignore", divide="ignore"):
out = (s / n).astype(np.float32)
out[~np.isfinite(arr)] = np.nan
return out
def describe(v):
v = np.asarray(v, dtype=float)
v = v[np.isfinite(v)]
if v.size == 0:
return None
return dict(n=v.size, mean=v.mean(), std=v.std(ddof=1) if v.size > 1 else np.nan,
median=np.median(v), p5=np.percentile(v, 5), p95=np.percentile(v, 95))
def fmt(name, d):
if d is None:
return f"{name}: nessun valore valido"
return (f"{name}: n={d['n']} media={d['mean']:.5f} std={d['std']:.5f} "
f"mediana={d['median']:.5f} p5={d['p5']:.5f} p95={d['p95']:.5f}")
def main():
ap = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("raster", help="GeoTIFF della profondita' di banda")
ap.add_argument("points", help="shapefile puntuale")
ap.add_argument("-o", "--output", help="file di uscita (.gpkg o .csv) con il valore per punto")
ap.add_argument("--radius", type=int, default=0,
help="raggio in pixel della finestra attorno al punto (0 = solo il pixel del punto)")
ap.add_argument("--radius-stat", choices=["max", "mean"], default="max",
help="statistica nella finestra se --radius > 0 (default: max)")
args = ap.parse_args()
with rasterio.open(args.raster) as ds:
arr = ds.read(1).astype(np.float32)
if ds.nodata is not None and np.isfinite(ds.nodata):
arr[arr == ds.nodata] = np.nan
arr[~np.isfinite(arr)] = np.nan
H, W = arr.shape
gdf = gpd.read_file(args.points)
n_total = len(gdf)
if gdf.crs is None:
sys.exit("[ERRORE] Lo shapefile non ha CRS (.prj mancante).")
gdf = gdf.explode(index_parts=False).reset_index(drop=True)
if not (gdf.geom_type == "Point").all():
sys.exit("[ERRORE] Lo shapefile deve contenere solo punti.")
if ds.crs is not None and gdf.crs != ds.crs:
print(f"Riproiezione dei punti da {gdf.crs} a {ds.crs}")
gdf = gdf.to_crs(ds.crs)
# coordinate -> indici di pixel con la trasformata inversa (indipendente
# dalla versione di rasterio, che con gli array a volte fallisce)
inv = ~ds.transform
x = gdf.geometry.x.values
y = gdf.geometry.y.values
cols = np.floor(inv.a * x + inv.b * y + inv.c).astype(int)
rows = np.floor(inv.d * x + inv.e * y + inv.f).astype(int)
inside = (rows >= 0) & (rows < H) & (cols >= 0) & (cols < W)
grid = neighborhood(arr, args.radius, args.radius_stat)
value = np.full(len(gdf), np.nan, dtype=np.float32)
value[inside] = grid[rows[inside], cols[inside]]
valid = np.isfinite(value)
print(f"Punti nello shapefile: {n_total} (dopo explode: {len(gdf)})")
print(f"Punti dentro l'immagine: {int(inside.sum())}")
print(f"Punti su pixel validi: {int(valid.sum())}")
if args.radius:
print(f"Finestra: {2 * args.radius + 1}x{2 * args.radius + 1} pixel, statistica = {args.radius_stat}")
print()
d_pts = describe(value[valid])
print(fmt("Per punto ", d_pts))
# punti che cadono nello stesso pixel: contarli una volta sola evita duplicati
key = rows[valid] * W + cols[valid]
_, first = np.unique(key, return_index=True)
d_pix = describe(value[valid][first])
print(fmt("Per pixel unico ", d_pix))
d_all = describe(grid[np.isfinite(grid)])
print(fmt("Tutti i pixel ", d_all))
if d_pix and d_all and d_all["std"] > 0:
eff = (d_pix["mean"] - d_all["mean"]) / d_all["std"]
print(f"\nScarto della media (pixel unici) rispetto a tutta l'immagine: "
f"{eff:+.2f} deviazioni standard dell'immagine")
if args.output:
out = gdf.copy()
out["depth"] = value
out["row"] = rows
out["col"] = cols
out["inside"] = inside
out["valid"] = valid
if args.output.lower().endswith(".csv"):
out.drop(columns="geometry").assign(x=gdf.geometry.x, y=gdf.geometry.y) \
.to_csv(args.output, index=False)
else:
out.to_file(args.output, driver="GPKG")
print(f"\nSalvato: {args.output}")
if __name__ == "__main__":
main()
Dai risultati si vede che i punti si confondono con il background
Il problema con questo approccio e' che il rumore domina il segnale oltre al fatto che l'assorbimento del crisotilo e' vicino a quello del carbonato del cemento
Nessun commento:
Posta un commento