In "Fundamental of geological and environmental remote sensing R.K.Vincent, Prentice Hall (1997)" cap 7 viene indicato che lo spostamento del Red Edge della vegetazione puo' essere un indicatore diagnostico per il tipo di mineralogia del suolo quando questo non possa essere visibile a causa della copertura vegetale
 |
| pag 197 del libro |
Questo spostamento e' indicativamente verso il blu per zone a maggiore mineralizzazione (anche se puo' accadere l'opposto)
 |
| Gavorrano |
Ho provato con lo script sottostante a verificare cosa accadeva nei dintorni dei depositi minerari della Toscana Meridionale
Questi sono gli spettri di alcuni punti nell'intorno
se si effettua la normalizzazione fondamentalmente non si vede nessun spostamento del red edg
plottando la posizione del red edge con come semplice spettro ma come mappa di ogni singolo pixel con lo script sottostante (per fare un lavoro completo avrei dovuto mascherare i punti che non sono vegetazione) si conferma che non ci sono particolari indicazion di variazioni del red edge attorno al sito minerario
import h5py
import numpy as np
import rasterio
from rasterio.transform import from_origin
file_path = 'gavorrano.h5'
output_tiff = 'gavorrano_REP_normalized.tif'
# Lunghezze d'onda nominali EnMAP VNIR (approssimative se mancano i metadati)
# In un file standard: Band 42 ~671nm, Band 52 ~703nm, Band 62 ~736nm, Band 75 ~779nm
def get_band_data(f, nr):
return f[f'bands/band_{nr:03d}'][:]
with h5py.File(file_path, 'r') as f:
band_indices = range(38, 81)
first_band = get_band_data(f, band_indices[0])
rows, cols = first_band.shape
cube = np.zeros((rows, cols, len(band_indices)), dtype=np.float32)
print("Caricamento e normalizzazione bande...")
for i, b_nr in enumerate(band_indices):
cube[:, :, i] = get_band_data(f, b_nr)
# 2. Normalizzazione Min-Max (0-1) lungo l'asse spettrale
c_min = cube.min(axis=2, keepdims=True)
c_max = cube.max(axis=2, keepdims=True)
denom_norm = c_max - c_min
# Evitiamo divisioni per zero
norm_cube = np.divide(cube - c_min, denom_norm,
out=np.zeros_like(cube),
where=denom_norm != 0)
r670 = norm_cube[:, :, 42-38] # Banda 42
r700 = norm_cube[:, :, 52-38] # Banda 52
r740 = norm_cube[:, :, 62-38] # Banda 62
r780 = norm_cube[:, :, 75-38] # Banda 75
re_reflectance = (r670 + r780) / 2
denom_rep = r740 - r700
rep = np.full((rows, cols), -999, dtype='float32')
mask = (denom_rep != 0)
rep[mask] = 700 + 40 * ((re_reflectance[mask] - r700[mask]) / denom_rep[mask])
rep = np.where((rep > 680) & (rep < 760), rep, -999)
west = 500000.0 # Min Easting
north = 5200000.0 # Max Northing
pixel_size = 30.0 # Resolution in meters
transform = from_origin(west, north, pixel_size, pixel_size)
#transform = from_origin(680000, 4750000, 30, 30)
crs = 'EPSG:32632'
# 6. Scrittura GeoTIFF
with rasterio.open(
output_tiff, 'w',
driver='GTiff',
height=rows, width=cols,
count=1, dtype='float32',
crs=crs,
transform=transform,
nodata=-999
) as dst:
dst.write(rep, 1)
print(f"Completato! File salvato: {output_tiff}")
