Visualizzazione post con etichetta Sentinel. Mostra tutti i post
Visualizzazione post con etichetta Sentinel. Mostra tutti i post

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


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)

domenica 14 aprile 2024

Tensorflow segmentazione UNET per frane

Ho trovato su Kaggle questo post in cui si applica la segmentazione tramite rete neurale UNET a corpi idrici

In estrema sintesi al contrario di Yolo la segmentazione non avviene tramite un box attorno all'oggetto di interesse ma tramite una maschera di pixel di forma arbitraria

Funziona per le frane? (ovviamente non sono il primo a pensarci..volevo solo provare). Ho trovata su Kaggle questo dataset in cui sono gia' disponibili le immagini e le maschere (e' un sottoinsieme di un dataset creato per la competizione Landslide4sense 2022 molto piu' voluminoso 3Gb che si trova a questo indirizzo)


Nel dataset ci sono due coppie di immagini. La prima e' un tassello Sentinel 2 truecolor di 128x128 pixel

La seconda una immagine maschera sempre di 128x128 pixel di tipo binario (frana/no frana) realizzata con intervento umano di fotointerpretazione come verita' a terra



La disposizione dei colori e' particolare come si vede dall'istogramma. Ogni colore e' una classe e le classi sono consecutive a partire dal nero. Se si apre il file maschera in Gimp non si vede niente se non si modifica la curva colori



ho applicato, usando Colab e la GPU allegata, lo script per i corpi idrici al dataset delle frane

L'unica modifica che e' ho fatto sul dataset e' stato rinominare le maschere. Il nome del file nella cartella Images deve essere identico a quello della corrispettiva cartella Masks



from functools import partial
import os

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tensorflow as tf
from tensorflow import keras




images_dir = 'water/Images'
masks_dir = 'water/Masks'

dirname, _, filenames = next(os.walk(images_dir))



@tf.function
def load_img_with_mask(image_path, images_dir: str = 'Images', masks_dir: str = 'Masks',images_extension: str = 'jpg', masks_extension: str = 'jpg') -> dict:
image = tf.io.read_file(image_path)
image = tf.image.decode_jpeg(image, channels=3)

mask_filename = tf.strings.regex_replace(image_path, images_dir, masks_dir)
mask_filename = tf.strings.regex_replace(mask_filename, images_extension, masks_extension)
mask = tf.io.read_file(mask_filename)
mask = tf.image.decode_image(mask, channels=1, expand_animations = False)
return (image, mask)

n_examples = 3
examples = [load_img_with_mask(os.path.join(images_dir, filenames[i])) for i in range(n_examples)]

fig, axs = plt.subplots(n_examples, 2, figsize=(14, n_examples*7), constrained_layout=True)
for ax, (image, mask) in zip(axs, examples):
ax[0].imshow(image)
ax[1].imshow(mask)
plt.show()



@tf.function
def resize_images(images, masks, max_image_size=1500):
shape = tf.shape(images)
scale = (tf.reduce_max(shape) // max_image_size) + 1
target_height, target_width = shape[-3] // scale, shape[-2] // scale
images = tf.cast(images, tf.float32)
masks = tf.cast(masks, tf.float32)
if scale != 1:
images = tf.image.resize(images, (target_height, target_width), method=tf.image.ResizeMethod.NEAREST_NEIGHBOR)
masks = tf.image.resize(masks, (target_height, target_width), method=tf.image.ResizeMethod.NEAREST_NEIGHBOR)
return (images, masks)

@tf.function
def scale_values(images, masks, mask_split_threshold = 128):
images = tf.math.divide(images, 255)
masks = tf.where(masks > mask_split_threshold, 1, 0)
return (images, masks)

@tf.function
def pad_images(images, masks, pad_mul=16, offset=0):
shape = tf.shape(images)
height, width = shape[-3], shape[-2]
target_height = height + tf.math.floormod(tf.math.negative(height), pad_mul)
target_width = width + tf.math.floormod(tf.math.negative(width), pad_mul)
images = tf.image.pad_to_bounding_box(images, offset, offset, target_height, target_width)
masks = tf.cast(tf.image.pad_to_bounding_box(masks, offset, offset, target_height, target_width), tf.uint8)
return (images, masks)

batch_size = 32
test_set_size = 300
validation_set_size = 250

dataset = tf.data.Dataset.list_files(images_dir + '/*.jpg', seed=42)

test_dataset = dataset.take(test_set_size)
dataset = dataset.skip(test_set_size)
test_dataset = test_dataset.map(load_img_with_mask)
test_dataset = test_dataset.map(scale_values)
test_dataset = test_dataset.shuffle(20)
test_dataset = test_dataset.map(lambda img, mask: resize_images(img, mask, max_image_size=2500))
test_dataset = test_dataset.map(pad_images)
test_dataset = test_dataset.batch(1).prefetch(5)


validation_dataset = dataset.take(validation_set_size)
train_dataset = dataset.skip(validation_set_size)
validation_dataset = validation_dataset.map(load_img_with_mask)
validation_dataset = validation_dataset.map(scale_values)
validation_dataset = validation_dataset.shuffle(20)
validation_dataset = validation_dataset.map(resize_images)
validation_dataset = validation_dataset.map(pad_images)
validation_dataset = validation_dataset.batch(1).prefetch(5)

train_dataset = train_dataset.map(load_img_with_mask)
train_dataset = train_dataset.map(scale_values)
train_dataset = train_dataset.shuffle(20)
train_dataset = train_dataset.map(resize_images)
train_dataset = train_dataset.map(pad_images)
train_dataset = train_dataset.batch(1).prefetch(5)

def get_unet(hidden_activation='relu', initializer='he_normal', output_activation='sigmoid'):
PartialConv = partial(keras.layers.Conv2D,
activation=hidden_activation,
kernel_initializer=initializer,
padding='same')
# Encoder
model_input = keras.layers.Input(shape=(None, None, 3))
enc_cov_1 = PartialConv(32, 3)(model_input)
enc_cov_1 = PartialConv(32, 3)(enc_cov_1)
enc_pool_1 = keras.layers.MaxPooling2D(pool_size=(2, 2))(enc_cov_1)
enc_cov_2 = PartialConv(64, 3)(enc_pool_1)
enc_cov_2 = PartialConv(64, 3)(enc_cov_2)
enc_pool_2 = keras.layers.MaxPooling2D(pool_size=(2, 2))(enc_cov_2)
enc_cov_3 = PartialConv(128, 3)(enc_pool_2)
enc_cov_3 = PartialConv(128, 3)(enc_cov_3)
enc_pool_3 = keras.layers.MaxPooling2D(pool_size=(2, 2))(enc_cov_3)
# Center
center_cov = PartialConv(256, 3)(enc_pool_3)
center_cov = PartialConv(256, 3)(center_cov)
# Decoder
upsampling1 = keras.layers.UpSampling2D(size=(2, 2))(center_cov)
dec_up_conv_1 = PartialConv(128, 2)(upsampling1)
dec_merged_1 = tf.keras.layers.Concatenate(axis=3)([enc_cov_3, dec_up_conv_1])
dec_conv_1 = PartialConv(128, 3)(dec_merged_1)
dec_conv_1 = PartialConv(128, 3)(dec_conv_1)
upsampling2 = keras.layers.UpSampling2D(size=(2, 2))(dec_conv_1)
dec_up_conv_2 = PartialConv(64, 2)(upsampling2)
dec_merged_2 = tf.keras.layers.Concatenate(axis=3)([enc_cov_2, dec_up_conv_2])
dec_conv_2 = PartialConv(64, 3)(dec_merged_2)
dec_conv_2 = PartialConv(64, 3)(dec_conv_2)
upsampling3 = keras.layers.UpSampling2D(size=(2, 2))(dec_conv_2)
dec_up_conv_3 = PartialConv(32, 2)(upsampling3)
dec_merged_3 = tf.keras.layers.Concatenate(axis=3)([enc_cov_1, dec_up_conv_3])
dec_conv_3 = PartialConv(32, 3)(dec_merged_3)
dec_conv_3 = PartialConv(32, 3)(dec_conv_3)
output = keras.layers.Conv2D(1, 1, activation=output_activation)(dec_conv_3)
return tf.keras.Model(inputs=model_input, outputs=output)



model = get_unet()

optimizer = tf.keras.optimizers.Nadam()
model.compile(loss='binary_crossentropy', optimizer=optimizer)

model.summary()



early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5)
lr_reduce = tf.keras.callbacks.ReduceLROnPlateau(monitor='val_loss', factor=0.3, patience=3, verbose=1)

epochs = 80
history = model.fit(train_dataset, validation_data=validation_dataset, epochs=epochs, callbacks=[early_stopping, lr_reduce])


n_examples = 10

fig, axs = plt.subplots(n_examples, 3, figsize=(14, n_examples*7), constrained_layout=True)
for ax, ele in zip(axs, test_dataset.take(n_examples)):
image, y_true = ele
prediction = model.predict(image)[0]
prediction = tf.where(prediction > 0.5, 255, 0)
ax[0].set_title('Original image')
ax[0].imshow(image[0])
ax[1].set_title('Original mask')
ax[1].imshow(y_true[0])
ax[2].set_title('Predicted area')
ax[2].imshow(prediction)

plt.show()

meanIoU = tf.keras.metrics.MeanIoU(num_classes=2)
for ele in test_dataset.take(test_set_size):
image, y_true = ele
prediction = model.predict(image)[0]
prediction = tf.where(prediction > 0.5, 1, 0)
meanIoU.update_state(y_true[0], prediction)
print(meanIoU.result().numpy())



Questi sono alcuni confronti della validazione con a sinistra l'immagine Sentinel, al centro la maschera fotointerpretata da utente umano ed a destra la maschera di predizione della rete neurale






la accuratezza e' molto buona (forse troppa ...avro' fatto qualche errore?) pari a 0.98




mercoledì 3 aprile 2024

Tensorflow e firme spettrali

Questa e' sola una prova di fattibilita', non ha interesse di uso pratico

Era un po' di tempo che volevo provare ad usare una rete neurale per classificare le firme spettrali da telerilevamento

Per semplicita' ho preso una immagine Sentinel 2 e selezionato 90 Pin sulla mappa (3 categorie: sedimenti fluviali, suolo nudo, vegetazione) per estrarre le firme tramite Optical/Spectrum View


 

C'e' da osservare che SNAP  estrae solo le bande Sentinel 2 a risoluzione a 10 m. Il file csv e' del tipo  (ho provato a fare il resampling per estrarre tutta la firma spettrale mostrata in grafico ma senza successo)

490,560,665,842,Classe
0.0737,0.0992,0.1278,0.2956,1
0.064,0.0882,0.1038,0.2976,1
0.0682,0.0964,0.1086,0.2792,1
0.0674,0.0956,0.1142,0.2832,1
0.0772,0.1046,0.1386,0.2642,1 

 

A questo punto ho modificato questo esempio di tensorflow basato sul dataset Iris sostituendo le feature della morfologia dei fiori con la risposta spettrale


import pandas as pd
import numpy as np
import tensorflow as tf
import seaborn as sns
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense,Dropout
from sklearn.preprocessing import LabelEncoder , StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, confusion_matrix

df = pd.read_csv('spettri.csv')
print(df.head())

le = LabelEncoder()
df['classe'] = le.fit_transform(df['classe'])

X = df.drop(columns=['classe'])
y = df['classe']
print(y.head())
print(y.value_counts())

X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.25,shuffle=True,random_state=7)

sc = StandardScaler()
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)

print(y_train)
y_train = tf.keras.utils.to_categorical(y_train,num_classes=3)

def get_models():
model = Sequential([
Dense(units=32,input_shape=(4,),activation='relu'),
Dense(units=32,activation='relu'),
Dropout(0.5),
Dense(units=3,activation='softmax')
])
return model

model = get_models()
model.compile(optimizer='Adam',loss='categorical_crossentropy',metrics=['accuracy'])
model.summary()

model.fit(X_train,y_train,epochs=100, verbose=2)

prediction = model.predict(X_test)
prediction = np.argmax(prediction,axis=-1)
acury = accuracy_score(y_test,prediction)
print(acury)
cm = confusion_matrix(y_test,prediction)
print(cm)


In fase di test il perggior riconoscimento e' dell' 86%

Questa e' la matrice di confusione del risultato peggiore
 

[[7 1 0]
 [0 5 2]
 [0 0 8]]



lunedì 24 aprile 2023

Satellite Derived Bathymetry Mapping

Una applicazione dello script a questo link applicato direttamente in EO Browser


Versilia

Booca L'Arno

Lago di massaciuccoli

Golfo di Spezia Foce del Magra

in linea di principio lo script puo' essere semplificato dato che la batimetria puo' essere stimata in modo relativa con il rapporto log10(B3/B4) di Sentinel 2 (attenzione...Rrs non e' la Surface Reflectance che si trova nel prodotto L2 di Sentinel ma e' la Remote Sensing Reflectance)




in ogni caso si vede che l'influenza del carico di sedimenti influenza la stima della batimetria (vedi immagine successiva. Batimetria ripresa da progetto Camp Italy)





venerdì 21 aprile 2023

Correlazione tra Sentinel 5P e stazioni a terra

In questo post, usando gli script di Earth Engine visti in precedenza, ho provato a mettere in correlazione il dato telerilevato di Sentinel 5P con alcune stazioni di monitoraggio aria Arpat (dati pubblici sul sito www.arpat.toscana.it) 


Parametro NO2
Vi e' un ottimo accordo tra i dati medi mensili delle concentrazioni di NO2 da satellite ed a terra

la curva di correlazione non e' incoraggiante



Parametro CH4
I dati da satellite mostrano il gradiente di aumento del metano che e' coerente con il valore di crescita mondiale 




Parametro S02
Per il parametro S02 non si possono fare correlazioni perche' quando le concentrazioni aumentano (nel periodo autunnale) il satellite non riporta dati
Stazione Arpat via Bassi Firenze

Mappa SO2 Piana Fiorentina Giugno 2019
Sentinel 5P





Parametro O3
Per l'ozono ho molti dubbi...comunque esporti i dati si vede che c'e' un offset di circa 3 mesi tra i dati telerilevati e quelli a terra

Parametro CO
Questo parametro e' critico per la correlazione con la stazione a terra perche' a Firenze la stazione di misura non e' in posizione di fondo ma a bordo di una viabilita' principale



Parametro PM10
Il parametro PM10 puo' essere derivato da UV Aerosol nel dataset di Sentinel 5P
Qui il  confronto delle medie mensili tra i dati satellitari ed una stazione a terra




Aerosol marzo 2022



qui un riferimento bibliografico della correlazione tra dati satellitari ed a terra


Geologi

  E so anche espatriare senza praticamente toccare strada asfaltata