Red Edge Position Enmap

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}")