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

Frana M.Javello e Sentinel 1

Qualche settimana fa un conoscente mi ha mandato um video di una corona di frana sul monte Javello legata all'evento piovoso del 6 novembre

Ho provato a vedere se riuscivo a vedere l'estensione della frana dalle immagini di Sentinel 2. Sono riuscito ad inviduare un aspetto interessante ovvero che sul versante era presente una zona disboscata che e' stata tagliata nel 2023 ed i corrispondenza di questa area 

Immagine dell'area disboscata

Prima del taglio forestale

dopo il taglio forestale

Come si puo' vedere dopo l'evento l'impluvio risulta essere interessato da trasporto solido (e' la striscia bianca che scende verso valla a circa meta' dell'immagine leggermente a sinistra)

Nonostante cio' in ottico della frana nemmeno l'ombra. Ho provato, cosa che non avevo mai provato, ad utilizzare Sentinel 1.

Per fare change detection con immagini radar si devono cercare immagini appartenenti alla stessa orbita. il primo tentativo e' fallito perche' il sensore illuminava il versante sbagliato del M.Ferrato e quindi non si mostravano modifiche al versante

Ho preso quindi la coppia 

S1A_IW_GRDH_1SDV_20231031T171457_20231031T171522_051012_062696_0499.SAFE
S1A_IW_GRDH_1SDV_20231112T171456_20231112T171521_051187_062C9B_8EBB.SAFE

Ciascuna immagine e' stata trattata con in SNAP con

Radar/Apply Orbit
Radar/Radiometric/S1 Remove Thermal Noise
Radar/Radiometric/Calibrate (attivando sia Sigma0 che Beta0)
Radar/Speckle Filtering/Single Product Speckle Filtering
Radar/Geometric/Terrain Correction/Range Doppler Terrain Correction
Raster/Subset (per inquadrare solo l'area di interesse)
Raster/Data Conversion/Linear to/fromDb
Radar/Coregistration (per passare da due prodotti ad un prodotto solo)

Per il change detection sono utilizzati diversi algoritmi (sia basati su dati SLC che GRD) come . Ho provato il metodo piu' semplice che corrisponde al fare il rapporto (non la differenza) tra le bande Beta0  prima e dopo l'evento

rapporto beta0

scontornando manualmente le aree

Il metodo non e' particolarmente diagnostico. L'algoritmo indicava una frana anche a monte di Schignano ma persone del luogo mi hanno detto che non si sono verificati movimenti franosi. Diciamo che puo' essere una indicazione ma non altro

in effetti sembra che Sentinel 1 riesca a verificare delle zone di potenziale frana dove Sentinel 2 non riusciva a vedere. Spero di riuscire a trovare il tempo per un sopralluogo di persona 

1) Strategies for landslide detection using open-access synthetic aperture radar backscatter change in Google Earth Engine Alexander L. Handwerger , Shannan Y. Jones , Pukar Amatya  Hannah R. Kerner  , Dalia B. Kirschbaum , and Mong-Han Huang https://doi.org/10.5194/nhess-2021-283

2) Exploring event landslide mapping using Sentinel-1 SAR backscatter products Michele Santangelo, Mauro Cardinali, Francesco Bucci, Federica Fiorucci, Alessandro Cesare Mondini  Geomorphology Volume 397, 15 January 2022, 108021

DS0152 e geofono

Volevo provare l'oscilloscopio DS0152 e visto che ho in prestito un geofono ho guardato il segnale in uscita. Francamente pensavo che il sensore facesse uscire poche decine di millivolt ma come si vede dalla foto si passa da qualche decina di Volts fino a qualche Volt (ho letto anche fino a 9V ma con sollecitazioni meccaniche che non sono ragionevoli per un geofono)

 Anche la lettura del periodo del geofono (15 Hz) risulta coerente con le caratteristiche di fabbrica

Estrazione di piani da pointcloud

 Avevo gia' provato a fare una estrazione dei piani su un affioramento roccioso con PCL qualche tempo fa

Questa volta ho provato ad utilizzare la libreria Open3D con l'algoritmo Ransac per l'estrazione di piani multipli da una scansione fatta con il tablet Tango 

Per l'algoritmo e' sufficiente definire il numero massimo di piani attesi, la distanza di soglia ed il numero di iterazioni

I vari piani oltre ad essere visualizzati vengono anche estratti come ply

import sys
import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt
input_path="./"
dataname="maiano.ply"
pcd = o3d.io.read_point_cloud(input_path+dataname)
o3d.visualization.draw_geometries([pcd])

#ransac piani multipli
segment_models={}
segments={}
max_plane_idx=3
rest=pcd
for i in range(max_plane_idx):
    colors = plt.get_cmap("tab20")(i)
    segment_models[i], inliers = rest.segment_plane(distance_threshold=0.1,ransac_n=3,num_iterations=5000)
    segments[i]=rest.select_by_index(inliers)
    segments[i].paint_uniform_color(list(colors[:3]))
    rest = rest.select_by_index(inliers, invert=True)
    print("pass",i,"/",max_plane_idx,"done.")
o3d.io.write_point_cloud("0.ply", segments[0])
o3d.visualization.draw_geometries([segments[0]] + [rest])
o3d.io.write_point_cloud("1.ply", segments[1])
o3d.visualization.draw_geometries([segments[1]] + [rest])
o3d.io.write_point_cloud("2.ply", segments[2])
o3d.visualization.draw_geometries([segments[2]] + [rest])

Esplosione Hunga-Tonga-Hunga-Ha'apai registrata in Toscana

Il giorno 15 gennaio alle ore UTC 04:27 e' avvenuta l'esplosione del vulcano Hunga-Tonga-Hunga-Ha'apai alle isole Tonga

Le immagini dal satellite sono impressionanti ma ancora piu' impressionante e' che circa 16 ore quasi dalla parte opposta del mondo e a circa 17000 Km di distanza ovvero in Toscana e' stata registrato un sensibile aumento della pressione atmosferica dovuto alla propagazione dell'onda di pressione. La velocita' con tutte le approssimazioni introdotte e' quella del suono 

Prendendo dati di stazioni meteo pubbliche si puo' vedere 

Stazione di Molin del Piano (Fi)

Una anomalia del sensore?

Osservatorio Astronomico Margherita Hack

Aulla

Carrara

Bordighera

Correlazione tra serie tempo

 Un semplice script che calcola la correlazione tra due serie tempo (due estensimetri) facendo scorrere una finestra mobilew    

Come si vede il versante si muove con costanza in modo uniforme (i due estensimetri hanno dati correlabili)...quando c'e' una accelerazione questa e' locale e fa decorrelare i due punti di misura. In seguito si reinstaura il movimento generale

I dati meteo sono sostanzialmente inutili perche' in inverno le precipitazioni sono unicamente nevose. Nel resto dell'anno le precipitazioni non sembrano avere una particolare stagionalita'

import datetime as dt

import pandas as pd

import numpy as np

from matplotlib import pyplot as plt

from sklearn.preprocessing import PolynomialFeatures

from scipy import signal

## tutti i dati sono ricampionamti su base giornaliera

## i dati sono campionati ogni 20 oppure ogni 15 minuti

## i dati radar sono campionati ogni 70 minuti

################### METEO ################################

## dal 24/7/2019 00;00 al 1/12/2021 23:59

parser = lambda date: dt.datetime.strptime(date, '%d/%m/%Y %H:%M')

meteo = pd.read_csv('meteo/meteo.csv', sep=';', parse_dates=['Data'], index_col=['Data'], infer_datetime_format=True)

#effettua il resampling sui 180 minuti

meteo_d = meteo.resample('180T').sum()

meteo_d.to_csv("meteo_180.csv")

print(meteo_d.shape)

#################### RADAR ###############################

parser = lambda date: dt.datetime.strptime(date, '%d/%m/%Y %H:%M:%S.%f')

# legge il csv e lo mette in un dataframe

dati_s1 = pd.read_csv('disp_raw_PTS_P2.csv', sep=',', parse_dates=['Tempo'], index_col=['Tempo'], date_parser=parser)

dati_s2 = pd.read_csv('disp_raw_PTS_P30.csv', sep=',', parse_dates=['Tempo'], index_col=['Tempo'], date_parser=parser)

# ricampiona in base giornaliera

dati_s1_d = dati_s1.resample('180T').mean()

dati_s2_d = dati_s2.resample('180T').mean()

#estrae solo i dati per i quali c'e' anche il dato meteo

#dati_s1_d_f = dati_s1_d[246:1107]

#dati_s2_d_f = dati_s1_d[246:1107]

#salva su file

dati_s1_d.to_csv("s1.csv")

dati_s2_d.to_csv("s2.csv")

# finestra mobile

mw = []

x1 = np.arange(862)

for i in range(1950,8846):

set1 = dati_s1_d[i-160:i]

set2 = dati_s2_d[i-160:i]

mw_v = np.corrcoef(set1.iloc[:,0],set2.iloc[:,0])

mw.append(mw_v[0,1])

plt.figure()

plt.subplot(311)

plt.plot(meteo_d.index,dati_s1_d[1950:8846], color='green', label='P2')

plt.plot(meteo_d.index,dati_s2_d[1950:8846], color='blue', label='P30')

plt.legend()

plt.xlabel('Tempo')

plt.ylabel('Raw Radar (mm)')

plt.title('P2 P30')

plt.subplot(312)

plt.plot(meteo_d.index,mw, color='red')

plt.xlabel('Tempo')

plt.ylabel('Indice corr.')

plt.subplot(313)

plt.bar(meteo_d.index,meteo_d['Pluv.[mm]'])

plt.xlabel('Tempo')

plt.ylabel('Prec. (mm)')

plt.show()

#### ricerca del lag

print("########lag")

correlation = signal.correlate(dati_s1_d-np.mean(dati_s1_d), dati_s2_d - np.mean(dati_s2_d), mode="full")

lags = signal.correlation_lags(len(dati_s1_d), len(dati_s2_d), mode="full")

lag = lags[np.argmax(correlation)]

print(lag)

Estensimetro ed FFT

Una per applicare un filtro ai dati dell'estensimetro bassato su analisi FFT (vedi qui e precedenti)

Dati originali con curva di interpolazione polinomiale

 Detrend sottraendo la polinomiale

Analisi FFT

Il primo picco ha un periodo di circa 24 giorni. I successivi sono relativi al ciclo giornaliero ed armoniche di ordine superiore

applicando un filtro passa basso si ottiene il grafico filtrato invertendo la FFT

import pandas as pd

import numpy as np

from matplotlib import pyplot as plt

from sklearn.preprocessing import PolynomialFeatures

from scipy.fft import fft, ifft

# legge il csv e lo mette in un dataframe

dati = pd.read_csv('dati_small_prg.csv', sep=',', parse_dates=['Data'], index_col=['Data'])

# taglia la prima parte dei dati

dati = dati[1:8500]

# resample

#print("**********************************************")

#dati.info()

#dati.resample('60min').mean()

#dati.info()

#print("**********************************************")

####### regressione

import scipy

def f(x, a, b, c):

return a*x**2 + b*x + c

#return a*np.log(x)+b

#return a*x+b

def residual(p, x, y):

return y - f(x, *p)

p0 = [1., 1., 1.]

x = np.arange(1,8500,1)

popt, pcov = scipy.optimize.leastsq(residual, p0, args=(x,dati['Est.[mm]']))

print("Parametri della regressione ")

print(popt)

print("Errore sui parametri della regressione")

print(pcov)

yn = f(x, *popt)

plt.plot(x, yn, 'or')

plt.plot(x, dati['Est.[mm]'])

plt.show()

#Calcolo dell'errore quadratico medio

#tra i dati veri e la curva di correlazione

# come radice quadrata della somma delle differenze al quadrato

# diviso per il numero di elementi

differenza = yn-dati['Est.[mm]']

err = np.square(yn-dati['Est.[mm]'])

s_err = np.sqrt(np.sum(err, axis=0)/err.shape)

print("Errore quadratico medio (mm)")

print(s_err)

np.savetxt('correla.csv', yn, delimiter=',', fmt='%f')

# cancella le colonne che non servono

dati.drop(dati.columns[[0,1,4,5]], axis=1, inplace=True)

print(dati.info())

#salva i dati sul csv

dati.to_csv('./nuovo.csv')

#https://docs.scipy.org/doc/scipy/reference/tutorial/fft.html

#https://pythonnumericalmethods.berkeley.edu/notebooks/chapter24.04-FFT-in-Python.html

plt.plot(differenza, label='Detrend')

plt.legend(loc='best')

plt.show()

#con il passo di campionamento di 1 ora c'e' una frequenza

# di 0.0008333 Hz

# corrispondente ad una 1 misura ogni 20 minuti 1/(60x20)

dati_fft = np.fft.fft(differenza)

N = len(dati_fft)

n = np.arange(N)

# sr sampling rate

# n numero di misure

sr = 0.00083333

T = N/sr

print(T)

freq = n/T

plt.plot(freq[1:4250], np.abs(dati_fft[1:4250]))

plt.xlabel('Freq (Hz)')

plt.ylabel('FFT Amplitude |X(freq)|')

plt.show()

#filtra

sig_fft_filtered = dati_fft.copy()

cut_off = 0.00083

sig_fft_filtered[np.abs(freq) < cut_off] = 0

filtered = ifft(sig_fft_filtered)

plt.plot(filtered)

plt.xlabel('Tempo')

plt.ylabel('mm')

plt.show()

################################################

# 4.87927E-7 (23.72 giorni)

# 1.16571E-5 (23.82 ore)

# 2.32256E-5 (2 armonica del secondo dato)

# 3.48999E-5 (3 armonica del secondo dato)

# 4.6367E-5 (4 armonica del secondo dato)

# 5.8022E-5 (5 armonica del secondo dato)

# 6.9479E-5 (6 armonica del secondo dato)

# 8.1181E-5 (7 armonica del secondo dato)

# 9.2730E-5 (8 armonica del secondo dato)

Autoencoder anomaly detection con tensorlfow

 Terzo tentativo di analisi dati di estensimetro con Tensorflow (iniziato qui). Si tratta di un adattamento dell'esempio sul sito di Keras

In questo post si cerca di impostare una anomaly detection mediante autoencoder

La prima anomalia nella serie dati e' in corrispondenza del movimento indicato dalla freccia nel grafico soprastante

I dati sono stati tagliati in modo da includere solo l'inizio dell'anomalia in modo da non istruire troppo la rete 

Il modello converge rapidamente con valori di loss e validation loss similari

sovrapponendo il modello ai dati di train si nota una ottima corrispondenza

sottrando i dati reali dal modello si possono estrapolare le anomalie. Indicato dalla freccia l'anomalia derivante dal movimento

di seguito il codice

# -*- coding: utf-8 -*-

"""timeseries_anomaly_detection_detrend3

Automatically generated by Colaboratory.

Original file is located at

https://colab.research.google.com/drive/12Kkjp_xazCmO4HrK0tmzVPyHIoYxJsYo

"""

import numpy as np

import pandas as pd

from tensorflow import keras

from tensorflow.keras import layers

from matplotlib import pyplot as plt

!rm detrend.*

!wget http://c1p81.altervista.org/detrend3.zip

!rm *.csv

!unzip detrend3.zip

df_small_noise=pd.read_csv(r'detrend3.csv', sep=':', header=0, low_memory=False, infer_datetime_format=True, parse_dates={'datetime':[0]}, index_col=['datetime'],usecols=['Data','detrend'])

print(df_small_noise.head())

print(df_small_noise.shape)

#df_small_noise = df_small_noise[:9500]

plt.plot(df_small_noise['detrend'])

plt.show()

# Normalize and save the mean and std we get,

# for normalizing test data.

training_mean = df_small_noise.mean()

training_std = df_small_noise.std()

df_training_value = (df_small_noise - training_mean) / training_std

print("Number of training samples:", len(df_training_value))

TIME_STEPS = 1000

# Generated training sequences for use in the model.

def create_sequences(values, time_steps=TIME_STEPS):

output = []

for i in range(len(values) - time_steps + 1):

output.append(values[i : (i + time_steps)])

return np.stack(output)

x_train = create_sequences(df_training_value.values)

print("Training input shape: ", x_train.shape)

model = keras.Sequential(

[

layers.Input(shape=(x_train.shape[1], x_train.shape[2])),

layers.Conv1D(

filters=32, kernel_size=7, padding="same", strides=2, activation="relu"

),

layers.Dropout(rate=0.2),

layers.Conv1D(

filters=16, kernel_size=7, padding="same", strides=2, activation="relu"

),

layers.Conv1DTranspose(

filters=16, kernel_size=7, padding="same", strides=2, activation="relu"

),

layers.Dropout(rate=0.2),

layers.Conv1DTranspose(

filters=32, kernel_size=7, padding="same", strides=2, activation="relu"

),

layers.Conv1DTranspose(filters=1, kernel_size=7, padding="same"),

]

)

model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001), loss="mse")

model.summary()

history = model.fit(

x_train,

x_train,

epochs=10,

batch_size=128,

validation_split=0.1,

callbacks=[

keras.callbacks.EarlyStopping(monitor="val_loss", patience=5, mode="min")

],

)

plt.plot(history.history["loss"], label="Training Loss")

plt.plot(history.history["val_loss"], label="Validation Loss")

plt.legend()

plt.show()

# Get train MAE loss.

x_train_pred = model.predict(x_train)

train_mae_loss = np.mean(np.abs(x_train_pred - x_train), axis=1)

plt.hist(train_mae_loss, bins=50)

plt.xlabel("Train MAE loss")

plt.ylabel("No of samples")

plt.show()

# Get reconstruction loss threshold.

threshold = np.max(train_mae_loss)

print("Reconstruction error threshold: ", threshold)

print(x_train.shape)

# Checking how the first sequence is learnt

plt.plot(x_train[288],label='Dati')

plt.plot(x_train_pred[288],label='Modello')

plt.legend()

plt.show()

anomalia = x_train[288] - x_train_pred[288]

plt.plot(anomalia)

plt.show()