Un ricordo dal passato (2010)
sabato 1 giugno 2024
venerdì 31 maggio 2024
Download and process Sentinel 5P data
Per scaricare in automatico i dati Sentinel ho utilizzato il progeto CDSETool
ho creato un file shape in coordinate WGS86 4236 nel folder dello script
from cdsetool.query import query_features, shape_to_wkt
from cdsetool.credentials import Credentials
from cdsetool.download import download_features
from cdsetool.monitor import StatusMonitor
from datetime import date
#chiavi di ricerca per Sentinel 5
# dict_keys(['maxRecords', 'index', 'page', 'identifier', 'geometry', 'box', 'lon', 'lat', 'radius', 'startDate', 'completionDate',
# 'productIdentifier', 'productType', 'processingLevel', 'platform', 'instrument', 'orbitNumber', 'sensorMode',
# 'updated', 'publishedAfter', 'publishedBefore', 'sortParam', 'sortOrder', 'status', 'exactCount', 'processingBaseline'])
PASSWORD_SENTINEL_ESA = "xxxxxxxx"
USER_SENTINEL_ESA = "l.innocenti@xxxxxx"
PATH_WHERE_TO_SAVE = "/home/luca/netcdf/"
geometry = shape_to_wkt("toscana.shp")
print(geometry)
features = query_features(
"Sentinel5P",
{
"startDate": "2024-05-30",
"completionDate": date(2024, 5, 31),
"maxRecords": 5,
#"processingLevel": "S5PL2",
"productType": "L2__O3____",
# title_ids : ["32TPP"], solo per Sentinel 2
"geometry": geometry
},
)
list(
download_features(
features,
PATH_WHERE_TO_SAVE,
{
"concurrency": 4,
"monitor": StatusMonitor(),
"credentials": Credentials(USER_SENTINEL_ESA, PASSWORD_SENTINEL_ESA),
},
)
)
ATTENZIONE: i files scaricati nonostante riportino una estensione nc in realta' sono zippati. Si devono quindi prima decomprimere e poi trattare come NetCDF.nc)
Per visualizzare in modo speditivo i dati si puo' usare Panoply. Alrimenti si puo' usare anche SNAP (per vedere i dati georiferiti https://www.youtube.com/watch?v=G8tVNbdu8-A)
il file nc puo' essere letto tramite xarray
in sintesi ci sono diversi array, uno contiene le misure, uno la longitudine del pixel ed uno la latitudine del pixel
La cosa piu' comoda per l'analisi e' salvare i dati in un database postgis. Attenzione che devono essere gestiti i valori nan che possono essere presenti
#!/usr/bin/env python
# coding: utf-8
import xarray as xr
import math
import psycopg
from shapely.geometry import LineString
from shapely import wkb
conn = psycopg.connect("dbname=postgis_db user=postgres password=Tvlgal55 host=localhost")
curs = conn.cursor()
#CREATE TABLE ozono(geom geometry, valore FLOAT not NULL, data TIMESTAMP without time zone not NULL)'
#data5p = xr.open_dataset('./S5P_NRTI_L2__O3_____20240529T122201_20240529T122701_34335_03_020601_20240529T130706/S5P_NRTI_L2__O3_____20240529T122201_20240529T122701_34335_03_020601_20240529T130706.nc', group='PRODUCT')
data5p = xr.open_dataset('./S5P_NRTI_L2__O3_____20240530T120201_20240530T120701_34349_03_020601_20240530T124545/S5P_NRTI_L2__O3_____20240530T120201_20240530T120701_34349_03_020601_20240530T124545.nc', group='PRODUCT')
#print(data5p)
#print(data5p.dims)
dimx= data5p.dims['scanline']
dimy=data5p.dims['ground_pixel']
tempo = data5p['ozone_total_vertical_column']['time']
data= "2024-05-30 12:02:01"
#print(dimx)
#print(dimy)
#print(tempo)
for x in range(0,dimx):
for y in range(0,dimy):
lon = str(data5p.longitude.values[0,x,y])
lat = str(data5p.latitude.values[0,x,y])
valore = data5p.ozone_total_vertical_column.values[0,x,y]
#print(str(data5p.longitude.values[0, x,y])+";"+str(data5p.latitude.values[0, x,y])+";"+str(data5p.ozone_total_vertical_column.values[0, x,y]))
if not(math.isnan(valore)):
curs.execute("INSERT INTO ozono (geom, valore,data) VALUES (ST_SetSRID(ST_MakePoint("+lon+","+lat+"), 4326),"+str(valore)+",'"+data+"')");
conn.commit()
SELECT avg(valore) as median FROM ozono WHERE ST_DistanceSphere(geom, ST_MakePoint(11.1,43.1)) <= 10000 GROUP BY data
Si puo' usare QGis per vestire i dati con un stile ed esportando il file .SLD
Il prodotto finale puo' essere distribuito con Openlayers su web
per indicizzare un campo spaziale in postgis si puo' usare la sintassi
i dati in geoserver possono essere filtrati tramite il campo CQL basandosi su un box oppure tramite un buffer attorno ad un punto
BBOX(geom, 10, 42.5, 12, 44.5)
DWithin(geom,POINT(11.38 43.3),140000,meters)
giovedì 23 maggio 2024
Regressione umidita' Random Forest con YDF
Volevo sperimentare la libreria YDF (Decision Trees) ed ho preso i dati del post precedente (soil moisture su terreni naturali)
Non ho fatto nessuna ottimizzazione, ho usato l'esempio base della regressione
(da notare che i grafici funzionano su Colab non riesco a visualizzarli su jupyter notebook locali)
nel codice sottostante e' stato attivato il GradientBoostTreeLearner ma vi sono anche RandomForestLearner, CartLearner
from google.colab import drive
drive.mount('/content/drive')
!pip install ydf
import ydf # Yggdrasil Decision Forests
import pandas as pd # We use Pandas to load small datasets
all_ds = pd.read_csv(f"/content/drive/MyDrive/soilmoisture.csv")
# Randomly split the dataset into a training (70%) and testing (30%) dataset
all_ds = all_ds.sample(frac=1)
split_idx = len(all_ds) * 7 // 10
train_ds = all_ds.iloc[:split_idx]
test_ds = all_ds.iloc[split_idx:]
# Print the first 5 training examples
train_ds.head(5)
model = ydf.GradientBoostedTreesLearner(label="soil_moisture",
task=ydf.Task.REGRESSION).train(train_ds)
evaluation = model.evaluate(test_ds)
print(evaluation)
evaluation
La sintassi e' estremamente concisa ma il risultato della regressione e' ottimale
per vedere quale ' la feature che influenza il modello si puo' usare
model.analyze(test_ds, sampling=0.1)
lunedì 20 maggio 2024
Regressione umidita' su dati iperspettrali con Tensorflow
Ho trovato su GitHub un dataset iperspettrale con data di taratura di umidita' su suoli naturali
Felix M. Riese and Sina Keller, “Introducing a Framework of Self-Organizing Maps for Regression of Soil Moisture with Hyperspectral Data,” in IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium, Valencia, Spain, 2018, pp. 6151-6154. (Link)
@inproceedings{riese2018introducing,
author = {Riese, Felix~M. and Keller, Sina},
title = {{Introducing a Framework of Self-Organizing Maps for Regression of Soil Moisture with Hyperspectral Data}},
booktitle = {IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium},
year = {2018},
month = {July},
address = {Valencia, Spain},
doi = {10.1109/IGARSS.2018.8517812},
ISSN = {2153-7003},
pages = {6151--6154},
}
Si tratta di circa 680 spettri con frequenze da 454 a 950 nm (125 bande). L'umidita' nel dataset varia da un minimo di 25.5% ad un massimo di 42.5% con un valore medio di 31.7%...purtroppo la variabilita' non e' enorme)
L'umidita' agisce sulla risposta spettrale non mediante un assorbimento localizzato su un picco ad una precisa lunghezza d'onda ma agisce su tutte le lunghezze d'onda abbasando la riflettanza all'aumentare dell'umidita'
Per questo motivo non e' banale fare una regressione tramite la rimozine del continuo
Ho provato ad usare lo script visto qui
In dottorato avevo visto che si otteneva una buona correlazione fissando una lunghezza d'onda e misurando la riflettanza dopo la rimozione del continuo ma come si vede dal grafico sottostante i coefficienti sono in funzione della lunghezza d'onda presa in considerazione
giovedì 16 maggio 2024
Unmixing iperspettrale di campioni di terreno naturale
Aggiornamento
Ho trovato che USGS Spectral library mette a disposizione anche degli spettri Fieldspec (2151 bande da 350 a 2500 nm) di campioni puri di specie mineralogiche e quindi sono stato in grado di avere gli spettri di tutti e 6 gli endmember mineralogici che descrivono per oltre il 95% la composizione dei campioni di suolo naturale. Ho provato a fare girare l'unmixing con tutte
Diciamo risultati piuttosto scarsi
-------------------------------------------------------------------------------------------
Sempre continuando a tirare fuori il materiale del dottorato da cui non avevo tirato fuori niente di significativo (vedi questo post e questo) ho ripreso la serie di spettri di campioni di terreno naturale del Mugello dei quali era stata fatta la caratterizzazione mineralogica quantitativa mediante XRD (analisi effettuata dal Dr. Lutterotti con il software https://luttero.github.io/maud/)
le specie mineralogiche erano indicate in 6 con a farla da padrone sono quarzo (26.7% wt medio) ed Illite (27.6% medio) seguiti da Montmorillonite (21.6% medio), Calcite (11.8% medio) e Kaolinite (3.07% medio)....da notare le concentrazioni perche' saranno determinanti per la discussione successiva
Dei campioni erano stati misurati spettri di laboratorio Fieldspec
Volevo vedere se era possibile fare l'unmixing spettrale mediante CLSU
Durante il dottorato avevo a disposizione dei campioni di standard naturale (illite, kaolinite, montmorillonite e calcite) ragionevolmente puri. Ho usato questi come endmember..il quarzo e' stato escluso perche' non ha assorbimenti significativi nel range del Fieldspec e l'albite e' stata esclusa per mancanza di standard (peraltro avrebbe anche assorbimenti molto modesti)
Per fare l'unmixing ho provato l'algoritmo Constrained LSU in SNAP. Con Octave mi sono create una immagine di 6x6 pixel e 2151 bande usando i primi 32 pixel per i campioni naturali e gli ultimi 4 per inserire di endmember per controllo
Dato che la somma di concentrazione in peso di calcite + montmorillonite + illite + kaolinite e' pari a circa il 43% (con un minimo del 29% ed un massimo del 60%) e visto che albite e quarzo non sono presi in considerazioni tra gli endmember dell'analisi spettrale sono state riscalate le concentrazioni di cal+mont+ill+kao a somma 100
I risultati non sono eccezionali..la correlazione della calcite e' similare a quanto ottenibile mediante la sola profondita' di picco a 2345 nm, per concentrazioni inferiori al 5% (come nel caso della Kaolinite) non c'e' nessun tipo di correlazione con il dato spettrale (questo aspetto e' stato uno dei risultati del dottorato)
lunedì 13 maggio 2024
domenica 12 maggio 2024
Unmixing spettri sintetici argille
AGGIORNAMENTO
ho provato ad usare questo metodo su campioni di miscele naturali di minerali argillosi (con concentrazioni determinate mediante XRD) ma ho miseramente fallito
Lo spettro dell'illite non viene estratto...anche forzando lo spettro puro nel dataset in modo da influenzare la VCA) non si ottengono risultati
=============================================
Per questo prova ho modificato lo script
https://github.com/ricardoborsoi/unmixing_spectral_variability
Non avendo Matlab gli script sono stati eseguiti in Octave con modestissime modifiche al codice per utilizzare i dati del fieldspec
Gli spettri sperimentali dei 3 endmember ssono stati ottenuti da campioni naturali acquistati con un grado di purezza mineralogica superiore al 90% (trattandosi di campioni naturali di minerali argillosi e' impossibile ottenere campioni puri) in condizioni di laboratorio
Da notare, nonostante siano tutti minerali argillosi, kaolinite e Montmorillonite hanno spettri simili mentre Illite e' molto dissimile.Si vedra' in seguito che cio' avra' ripercusssioni sull'unmixing
Esempio di spettro di sintesi
Lo script genera 2500 spettri di sintesi con differenti valori di concentrazione calcolate come modello Gaussiano random (mappa 50x50 pixel del 2151 bande)
La distribuzione delle concentrazioni dei 3 endmember e'risultata cosi' distribuita
Sugli spettri di sintesi viene applicato VCA per determinare gli endmembers
Gli spettri sono quindi elaborati mediante differenti algoritmi di unmixing (nella matrice A sono contenute le concentrazioni reali). Piu' punti stanno sulla retta y=x migliore e' l'algoritmo
ELMM
FCLS
MESMA
RUSAL
In conclusione l'unico algoritmo che e' riuscito ad effettuare l'unmixing corretto e' stato MESMA
Iscriviti a:
Post (Atom)
-
Per lo shield WiFly e' stata effettuata la patch descritta qui Per la connessione del lettore Rfid si usano i seguenti pin VCC ... -
Questo post e' a seguito di quanto gia' visto nella precedente prova Lo scopo e' sempre il solito: creare un sistema che permet... -
In questo post viene indicato come creare uno scatterplot dinamico basato da dati ripresi da un file csv (nel dettaglio il file csv e' c...