Primo spettro di riflettanza con SparkFun Triad Spectroscopy Sensor - AS7265x

 Lo spettro di riflettanza e' stato ottenuto dividendo il valore calibrato della foglia per il valore calibrato nel corrispondente canale dello spectralon

 

 

	spettro	spectralon	foglia
410	0.104046939572376	7631.94	794.08
435	0.100663330010594	1916.09	192.88
460	0.0941060471276098	5496.99	517.3
485	0.0846332588803228	2220.64	187.94
510	0.127077186925102	2647.21	336.4
535	0.209306506039295	3170.9	663.69
560	0.204982889816192	1385.14	283.93
585	0.171727663387711	1473.03	252.96
610	0.0880803424432005	4555.5	401.25
645	0.121946929339632	795.92	97.06
680	0.123669783940664	1054.34	130.39
705	0.278338126099833	232.99	64.85
730	0.766924924721724	295.57	226.68
760	0.718044428259286	219.68	157.74
810	0.575862515972518	539.99	310.96
860	1.02470169644788	2543.55	2606.38
900	0.777308820373417	139.79	108.66
940	0.344453330286759	87.53	30.15

Primi passi con SparkFun Triad Spectroscopy Sensor - AS7265x

 Primi approcci per vedere se il sensore multispettrale e' affidabile per l'uso multispettrale 

In questo prova e' stata verificata la ripetibilita' nel tempo delle misure ponendo il sensore su uno spectralon ed effettuando 350 misure di bianco (il supporto del sensore e' stato ripreso da qui https://github.com/Scottapotamas/AS7265x-triad-ui/tree/master/mechanical)

 

 

Lo script per la misura ha previsto il massimo tempo di integrazione (con valori bassi in IR praticamente non c'era segnale)

sensor.setGain(AS7265X_GAIN_16X);
sensor.setMeasurementMode(AS7265X_MEASUREMENT_MODE_6CHAN_ONE_SHOT);

sensor.setIntegrationCycles(255);

L'errore per ogni canale e' stato stimato prendendo la standard deviation e dividendo per la media

A	0.44%
B	0.49%
C	0.25%
D	0.34%
E	0.05%
F	0.05%
G	0.07%
H	0.07%
R	0.10%
I	0.07%
S	0.09%
J	0.08%
T	0.13%
U	0.00%
V	0.37%
W	0.74%
K	1.55%
L	0.56%

Come si deve l'errore e' piuttosto vario tra i vari canali

 

Plottando la misura di ogni canale nel tempo si osserva una chiara deriva per ogni canale con in alcuni casi l'accenno di un valore asintotivo mentre in altri no 

 

 

 

 

La Dark Current e' stabile (prova su 450 misure)

	A	B	C	D	E	F	G	H	R	I	S	J	T	U	V	W	K	L
St Dev	0.00	0.00	0.00	0.00	0.00	0.00	0.03	0.00	0.00	0.00	0.46	0.00	0.40	0.40	0.27	0.00	0.00	0.00
Media	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.41	1.02	0.00	0.80	0.83	0.84	1.08	0.64	0.00

Geologi

 

E so anche espatriare senza praticamente toccare strada asfaltata

Campagna verita' a terra Fieldspec

 Un ricordo dal passato (2010)

 

 

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()

Per effettuare una query geografica si puo' usare il motore di postgis indicando il punto geografico desiderato ed un buffer..con questa query i dati dello stesso giorno vengono mediati

 

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

 

La vestizione puo' essere importata in Geoserver

Il prodotto finale puo' essere distribuito con Openlayers su web

 

per indicizzare un campo spaziale in postgis si puo' usare la sintassi

CREATE INDEX idx_ozono ON ozono USING gist (geom);

 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)

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)

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  

 

dallo scatterplot dei dati realmente misurati e quelli del modello derivanti dagli spettri nel dataset di test si vede una buona correlazione
 

 

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