giovedì 12 febbraio 2026

Anomalie in NDVI con algoritmo Continuous Change Detection (CCDC)

Attenzione: questo plugin non funziona su Debian 13 per delle dipendenze non risolte. Funziona in modo corretto su  Debian 12 con QGis 3.22 LTS

Attenzione 2 : la funzione CCDC e' presente su Google Earth Engine solo sui profili a pagamento...per questo motivo ho usato la versione di QGis...Aggiornamento: a quanto sembra non e' necessario il profilo a pagamento ma che lo script sia collegato ad un progetto di Google Cloud ed autorizzato (un progetto autorizzato si riconosce dall'icona in alto a destra)

Progetto autorizzato

Progetto non autorizzato

 

 Questo plugin usa l'algorimo Continous Change Detection che si basa su una funzione armonica per evidenziare brusche variazioni del parametro analizzato




Questo e' lo script per Google Earth Engine

// 1. DEFINIZIONE ROI

var areaOfInterest = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.896754172154326, 43.02122523564496],
          [10.896754172154326, 43.002271538414895],
          [10.935463812657256, 43.002271538414895],
          [10.935463812657256, 43.02122523564496]]], null, false);


Map.centerObject(areaOfInterest, 15);

// 2. PREPARAZIONE DATI
var prepareData = function(image) {
  var optical = image.multiply(0.0000275).add(-0.2);
  var ndvi = optical.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
  return image.addBands(ndvi).select('NDVI');
};

var collection = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2")
  .merge(ee.ImageCollection("LANDSAT/LT05/C02/T1_L2"))
  .merge(ee.ImageCollection("LANDSAT/LE07/C02/T1_L2"))
  .filterBounds(areaOfInterest)
  .map(prepareData)
  .sort('system:time_start');

// 3. ESECUZIONE CCDC CON GESTIONE ERRORI
// Proviamo a chiamare l'algoritmo nativo in modo più "protetto"
var ccdc;
try {
  ccdc = ee.Algorithms.TemporalSegmentation.ccdc({
    collection: collection,
    breakThreshold: 1.0,
    minObservationsForSegment: 6
  });
  print('CCDC eseguito con successo:', ccdc);
} catch (e) {
  print('Errore critico: Il tuo account non ha accesso a CCDC tramite API diretta.', e);
  print('Soluzione: Assicurati che lo script sia collegato a un Cloud Project (in alto a destra).');
}

// 4. SE CCDC HA FUNZIONATO, VISUALIZZA
if (ccdc) {
  var magnitude = ccdc.select('NDVI_magnitude').arrayGet([0]);
  var tBreak = ccdc.select('tBreak').arrayGet([0]);

  var magViz = {min: -0.5, max: 0.5, palette: ['#d73027', '#f4a582', '#f7f7f7', '#d9f0d3', '#1a9850']};
  var dateViz = {min: 1984, max: 2024, palette: ['#0000ff', '#00ffff', '#ffff00', '#ff0000']};

  Map.addLayer(magnitude.clip(areaOfInterest), magViz, 'Magnitudo Cambiamento');
  Map.addLayer(tBreak.clip(areaOfInterest), dateViz, 'Anno del Cambiamento');
}

Anomalie in NDVI

Questo script calcola l'anomalia dell'NDVI rispetto ad un valore medio

Unica accortezza.: non sono stati presi dati di Landsat 8 e 9 perche' si fondono male con Landsat 7






Map.centerObject(areaOfInterest, 15);
// 2. Funzione per NDVI scalato
var getNDVI = function(image) {
  var optical = image.multiply(0.0000275).add(-0.2);
  return optical.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI')
    .copyProperties(image, ['system:time_start']);
};
// 3. Collezioni Landsat 4, 5, 7
var collection = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
  .merge(ee.ImageCollection("LANDSAT/LE07/C02/T1_L2"))
  .filterBounds(areaOfInterest)
  .map(getNDVI);
// 4. CREAZIONE DELLE DUE EPOCHE (Baseline vs Recente)
// Epoca 1: Il passato (es. 1984-1994)
var past = collection.filterDate('1984-01-01', '1994-12-31').median();
// Epoca 2: Il presente (es. 2012-2022)
var present = collection.filterDate('2012-01-01', '2022-12-31').median();
// 5. CALCOLO DELLA MAGNITUDO DEL CAMBIAMENTO
// Differenza tra oggi e ieri
var change = present.subtract(past);
// 6. VISUALIZZAZIONE
var changeViz = {
  min: -0.3,
  max: 0.3,
  palette: ['#d73027', '#f4a582', '#f7f7f7', '#d9f0d3', '#1a9850']
};
Map.addLayer(past.clip(areaOfInterest), {min:0, max:0.8}, 'NDVI Passato (Anni 80/90)', false);
Map.addLayer(present.clip(areaOfInterest), {min:0, max:0.8}, 'NDVI Recente (Anni 2000/10)', false);
Map.addLayer(change.clip(areaOfInterest), changeViz, 'Cambiamento Netto (Magnitudo)');
// 7. Esportazione PNG del Cambiamento
var exportImg = change.visualize(changeViz);
print('Link PNG Cambiamento Storico:', exportImg.getThumbURL({
  dimensions: 1024,
  region: areaOfInterest,
  format: 'png'
}));


mercoledì 11 febbraio 2026

Mappa NDVI medio annuale con Google Earth Engine

Al contrario del precedente post https://debiaonoldcomputers.blogspot.com/2026/02/analisi-multitemporale-ndvi-con-google.html qui l'idea e' di creare una mappa di NDVI al posto della serie tempo per vedere anomalie spaziali

 


 

Landsat 8 

 // 1. Define the Study Area

var areaOfInterest = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.904865172222156, 43.01615135702198],
          [10.904865172222156, 43.00648615894694],
          [10.925807860210437, 43.00648615894694],
          [10.925807860210437, 43.01615135702198]]], null, false);

// 2. Pre-processing Function for Landsat 8
var processL8 = function(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var ndvi = opticalBands.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
  return image.addBands(ndvi).select('NDVI');
};

// 3. Load 2025 Median NDVI
var ndvi2025 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterBounds(areaOfInterest)
  .filterDate('2025-01-01', '2025-12-31')
  .map(processL8)
  .median()
  .clip(areaOfInterest);

// 4. Visualization Parameters
var palette = [
  '#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#056201', '#011301'
];
var ndviViz = {min: 0, max: 0.9, palette: palette};

// 5. Display the Map
Map.addLayer(ndvi2025, ndviViz, 'NDVI Median 2025');

// --- 6. ADDING THE LEGEND ---

// Create the legend panel
var legend = ui.Panel({
  style: {
    position: 'bottom-left',
    padding: '8px 15px'
  }
});

// Create legend title
var legendTitle = ui.Label({
  value: 'NDVI Index (2025)',
  style: {fontWeight: 'bold', fontSize: '16px', margin: '0 0 4px 0', padding: '0'}
});
legend.add(legendTitle);

// Helper function to create rows in the legend
var makeRow = function(color, name) {
  var colorBox = ui.Label({
    style: {
      backgroundColor: color,
      padding: '8px',
      margin: '0 0 4px 0'
    }
  });
  var description = ui.Label({
    value: name,
    style: {margin: '0 0 4px 6px'}
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')
  });
};

// Define labels for the palette colors
var names = ['Low/No Veg (< 0.1)', 'Soil / Sparse', 'Moderate Veg', 'Dense Veg', 'Healthy Forest', 'Peak Vigour (> 0.8)'];
var colors = ['#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#011301'];

// Add rows to legend
for (var i = 0; i < 6; i++) {
  legend.add(makeRow(colors[i], names[i]));
}

// Add legend to the map
Map.add(legend);

Export.image.toDrive({
  image: ndvi2025,            // The variable name of your NDVI map
  description: 'NDVI_2025_Tuscany', 
  folder: 'EarthEngine_Exports', // Name of the folder in your Google Drive
  fileNamePrefix: 'NDVI_2025_Site',
  region: areaOfInterest,     // Export only the area inside your polygon
  scale: 30,                  // Landsat resolution is 30 meters
  crs: 'EPSG:4326',           // Standard Coordinate Reference System (WGS84)
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true      // Makes the file easier to open in QGIS/ArcGIS
  }
});

 

Landsat 7

// 1. Define the Study Area
var areaOfInterest = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.904865172222156, 43.01615135702198],
          [10.904865172222156, 43.00648615894694],
          [10.925807860210437, 43.00648615894694],
          [10.925807860210437, 43.01615135702198]]], null, false);
// 2. Pre-processing Function for Landsat 8
var processL8 = function(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var ndvi = opticalBands.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
  return image.addBands(ndvi).select('NDVI');
};

// 3. Load 2025 Median NDVI
var l7_2005 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
  .filterBounds(areaOfInterest)
  .filterDate('2005-01-01', '2005-12-31')
  .map(function(image) {
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
    // Landsat 7 NDVI: (B4 - B3) / (B4 + B3)
    var ndvi = opticalBands.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
    return image.addBands(ndvi).select('NDVI');
  })
  .median() // Create a clean composite for the year
  .clip(areaOfInterest);

// 4. Visualization Parameters
var palette = [
  '#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#056201', '#011301'
];
var ndviViz = {min: 0, max: 0.9, palette: palette};

// 5. Display the Map
Map.addLayer(l7_2005, ndviViz, 'NDVI Median 2005');

// --- 6. ADDING THE LEGEND ---

// Create the legend panel
var legend = ui.Panel({
  style: {
    position: 'bottom-left',
    padding: '8px 15px'
  }
});

// Create legend title
var legendTitle = ui.Label({
  value: 'NDVI Index (2025)',
  style: {fontWeight: 'bold', fontSize: '16px', margin: '0 0 4px 0', padding: '0'}
});
legend.add(legendTitle);

// Helper function to create rows in the legend
var makeRow = function(color, name) {
  var colorBox = ui.Label({
    style: {
      backgroundColor: color,
      padding: '8px',
      margin: '0 0 4px 0'
    }
  });
  var description = ui.Label({
    value: name,
    style: {margin: '0 0 4px 6px'}
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')
  });
};

// Define labels for the palette colors
var names = ['Low/No Veg (< 0.1)', 'Soil / Sparse', 'Moderate Veg', 'Dense Veg', 'Healthy Forest', 'Peak Vigour (> 0.8)'];
var colors = ['#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#011301'];

// Add rows to legend
for (var i = 0; i < 6; i++) {
  legend.add(makeRow(colors[i], names[i]));
}

// Add legend to the map
Map.add(legend);

Export.image.toDrive({
  image: l7_2005,                // Updated to the 2005 variable
  description: 'NDVI_L7_2005_Tuscany', 
  folder: 'EarthEngine_Exports', 
  fileNamePrefix: 'NDVI_Landsat7_2005',
  region: areaOfInterest,     
  scale: 30,                  
  crs: 'EPSG:4326',           
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true      
  }
});


 

 

Analisi multitemporale NDVI con Google Earth Engine

Qualche esperimento con Google Earth Engine per vedere come si modifica nel tempo NDVI (arco di decenni fornito da Landsat)

 


 

// 1. Coordinate e Geometria
//var lon = 10.913507726271437;
//var lat = 43.01194709215584;

var lon = 10.926468160041452;
var lat = 43.01919548868449;


var poi = ee.Geometry.Point([lon, lat]);

// 2. Funzione di pre-elaborazione (Mask Cloud e Scaling)
// Valida per Landsat 5 e Landsat 7 (Collection 2 Level 2)
function preprocessLandsat(image) {
  var qa = image.select('QA_PIXEL');
  
  // Maschera per nuvole, ombre e cirri
  var mask = qa.bitwiseAnd(1 << 1).eq(0)
    .and(qa.bitwiseAnd(1 << 2).eq(0))
    .and(qa.bitwiseAnd(1 << 3).eq(0))
    .and(qa.bitwiseAnd(1 << 4).eq(0));

  // Fattori di scala per Collection 2
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  
  // Calcolo NDVI (NIR è B4, RED è B3 per entrambi i sensori L5 e L7)
  var ndvi = opticalBands.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
  
  return image.addBands(ndvi).updateMask(mask);
}

// 3. Caricamento Collezione LANDSAT 5 (1985 - 2010)
var l5Col = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
  .filterBounds(poi)
  .filterDate('1985-01-01', '2010-01-01')
  .map(preprocessLandsat);

// 4. Caricamento Collezione LANDSAT 7 (2010 - 2024)
var l7Col = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
  .filterBounds(poi)
  .filterDate('2010-01-01', '2024-01-01')
  .map(preprocessLandsat);

// 5. Unione delle due collezioni
var mergedCol = l5Col.merge(l7Col);

// 6. Creazione del Grafico Temporale
var chart = ui.Chart.image.series({
  imageCollection: mergedCol.select('NDVI'),
  region: poi,
  reducer: ee.Reducer.mean(),
  scale: 30,
  xProperty: 'system:time_start'
})
.setOptions({
  title: 'Serie Storica NDVI: Landsat 5 e Landsat 7 (1985-2024)',
  vAxis: {title: 'NDVI', viewWindow: {min: -0.1, max: 1}},
  hAxis: {title: 'Anno', format: 'YYYY'},
  lineWidth: 1,
  pointSize: 2,
  trendlines: {0: {color: 'red', lineWidth: 1, opacity: 0.7}}, // Linea di tendenza globale
  series: {
    0: {color: '228B22'} // Colore verde foresta
  }
});

// 7. Visualizzazione
print(chart);
Map.centerObject(poi, 13);
Map.addLayer(poi, {color: 'red'}, 'Sito di analisi'); 

 

 

uno script piu' evoluto si puo'trovare qui. L'interesse e' che non viene calcolato un singolo punto ma il valore medio di un'area che risulta sicuramente piu' significativo
 

modificando per l'area di interesse

 

// 1. Definizione dell'Area di Studio (Rettangolo)
var areaOfInterest = ee.Geometry.Polygon([
  [[10.910535746643099, 43.01138938313477],
   [10.91448395831302, 43.01138938313477],
   [10.91448395831302, 43.00825125524326],
   [10.910535746643099, 43.00825125524326],
   [10.910535746643099, 43.01138938313477]]
]);

Map.centerObject(areaOfInterest, 15);
Map.addLayer(areaOfInterest, {color: 'd63000'}, 'Area di Studio');

// 2. Funzione Processamento con Soglia NDVI > 0.55
var processLandsat = function(image, nir, red) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var ndvi = opticalBands.normalizedDifference([nir, red]).rename('NDVI');
  
  // FILTRO: Manteniamo solo i pixel dove NDVI è maggiore di 0.55
  var maskThreshold = ndvi.gt(0.55);
  
  return image.addBands(ndvi).select('NDVI')
    .updateMask(maskThreshold) // Rimuove i valori bassi
    .copyProperties(image, ['system:time_start']);
};

// 3. Caricamento Collezioni
var L5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B4', 'SR_B3'); });

var L7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B4', 'SR_B3'); });

var L8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B5', 'SR_B4'); });

var L9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B5', 'SR_B4'); });

var mergedNDVI = L5.merge(L7).merge(L8).merge(L9);

// 4. Composizione Mensile (Massimo Valore)
var years = ee.List.sequence(1984, 2024);
var months = ee.List.sequence(1, 12);

var monthlyNDVI = ee.ImageCollection.fromImages(
  years.map(function(y) {
    return months.map(function(m) {
      var d = ee.Date.fromYMD(y, m, 1);
      return mergedNDVI
        .filter(ee.Filter.calendarRange(y, y, 'year'))
        .filter(ee.Filter.calendarRange(m, m, 'month'))
        .max() // Seleziona il valore più alto (più verde) del mese
        .set('system:time_start', d.millis());
    });
  }).flatten()
);

// Rimuoviamo i periodi senza dati (es. inverno o periodi sotto soglia)
var monthlyNDVI_Clean = monthlyNDVI.filter(ee.Filter.listContains('system:band_names', 'NDVI'));

// 5. Grafico
var chart = ui.Chart.image.series({
  imageCollection: monthlyNDVI_Clean,
  region: areaOfInterest,
  reducer: ee.Reducer.mean(),
  scale: 30,
  xProperty: 'system:time_start'
})
.setOptions({
  title: 'NDVI Filtrato (>0.55) - Peak Vegetation Trend',
  vAxis: {title: 'NDVI', viewWindow: {min: 0.55, max: 0.9}},
  hAxis: {title: 'Anno', format: 'YYYY'},
  trendlines: {0: {color: 'red', lineWidth: 1}},
  series: {0: {color: '228B22', pointSize: 2, lineWidth: 1}}
});

print(chart); 

 


 


in generale sembra che NDVI sia migliorato a partire dal 2013

domenica 8 febbraio 2026

Flauto traverso barocco stampato 3D

Ho trovato questo modello di fluato traverso barocco da stampare https://www.thingiverse.com/thing:4572430 che e' una copia di un modello del 1740 di Bizet

  


I tenoni sono abbastanza sprecisi...ero partito per usare il nastro di teflon per garantire la tenuta ma alla fine per fare spessore ho usato in nastro carta. Per il tampone della chiave del Reb ho usato nastro biadesivo 

 

 



Io non saro' certo un grande flautista ma io non riesco a suonarlo. Non c'e' volume 

Inoltre la distanza dei fori e' diabolicamente molto superiore a quella di un flauto moderno. Peraltro la diteggiatura barocca, spesso con forchette, non aiuta 


 Esiste un altro modello disponibile da stampare che e' una copia di un flauto del 1750 conservato ad Amsterdam ed opera di Deppe  https://www.rijksmuseum.nl/en/collection/object/Flute--6e8161fcef916a7207f25a0812507870

Di questo secondo flauto esiste un video con nei commenti il link per i file STL.  

 


 

 

 

lunedì 2 febbraio 2026

Correzione atmosferica su dati EnMap

 Mi e' stato girato un articolo  

Detecting rare earth elements using EnMAP hyperspectral satellite data: a case study from Mountain Pass, California

in cui si indica che e' possibile mappare il neodimio da satellite usando bande VNIR. Mi sono incuriosito ed ho provato a scaricare i dati pubblici di Enmap relativi alla piu' grande miniera di Bayan Obo in Cina

Longitude:    109°54'55" E    degree
Latitude:    41°40'59" N    degree

 Dal punto di vista degli assorbimenti si vede che gli assorbimenti significativi sono a circa 740, 798 ed 870 nm (per rimanere nel campo del VNIR)

Ho provato a selezionare alcuni punti sia all'interno dell'area di escavazione che all'esterno
 
 

e questi sono gli spettri


 come si vede c'e' un picco di assorbimento a 760 nm ma questo e' presente sia nei dati dentro alla miniera che all'esterno ed i picchi sono via via che lo spettro ha mediamente una riflettanza maggiore

piu' in dettaglio 



 come si vede nella zona di 760 c'e' un profondo picco di assorbimento dell'ossigeno.

I dati che ho scaricato da Enmap erano di livello L2A quindi gia' trattati per la rimozione del contributo atmosferico  

L'ipotesi che l'algoritmo di rimozione atmosferica non abbia eliminato completamente il contributo dell'ossigeno e che quindi i dati non siano significativi nell'intorno della lunghezza d'onda di 760 nm 

 

 

 


domenica 1 febbraio 2026

Covox Speech Thing

Nel cassonetto delle pile usate a lavoro ho trovato questo circuito ed e' stato un tutto negli anni 80

Non lo ho mai usato (sono partito da un clone Soundblaster ad 8 bit) ma questa e' una primordiale scheda audio da collegare alla porta parallela

  


Originariamente nota come Covox Speech Thing (questa e' un clone) e' di fatto un convertitore Digitale Analogico ad 8 Bit basato solo resistenze in configurazione R-2R (https://en.wikipedia.org/wiki/Resistor_ladder)

Al termine c'e' uno stadio di amplificazione (da qui la necessita' delle pile) e l'uscita audio..nessun integrato..talmente semplice ed economico che se non ricordo male era incluso nelle confezioni di alcuni giochi

 Il problema di questi aggeggi e' che non avendo un timer se il processore del computer era impiegato in altri compiti oltre a quelli di aggiornare il suono la musica rallentava 

 

 

 



 

 

 

 

 

 

Anomalie in NDVI con algoritmo Continuous Change Detection (CCDC)

Attenzione: questo plugin non funziona su Debian 13 per delle dipendenze non risolte. Funziona in modo corretto su  Debian 12 con QGis 3.22 ...