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

martedì 19 novembre 2019

Riconoscimento foraminiferi con rete neurale Tensorflow Lite

Questo e' sempre stato uno dei miei sogni quando studiavo micropaleontologia....avere un sistema automatico (o semi automatico) che mi aiutasse nel riconoscimento di foraminiferi (planctonici del Miocene nel caso specifico)



Con Tensorflow Lite si puo' arrivare a qualcosa di simile
Prima di iniziare pero e' necessario avere una buona base di immagini....in cio' viene in aiuto il sito
http://www.endlessforams.org/ da cui e' possibile scaricare le fotografie di foraminiferi gia' classificati

Nello specifico, per il mio test, ho selezionato le prime 100 foto delle specie

Globigerina Falconensis
Globigerina Bulloides
Globigerinella Calida
Globigerinella Siphonifera
Globorotalia Crassiformis
Orbulina Universa

alcune sono state scelte volutamente con morfologia simile in modo da vedere il grado di risoluzione della macchina neurale. Queste 600 immagini sono in train dataset mentre sono state selezionate ulteriore 12 immagini (2 per specie) come test dataset; le immagini di test ovviamente non sono incluse nel set di apprendimento

Le immagini sono organizzate in una struttura di directory come in immagine seguente

In pratica il nome della directory e' la label utilizzata da Tensorflow per image recognition

E' giunto il momento di creare la libreria personalizzata di Tensorflow con il seguente script Python

Si apre il virtual environment con venv/bin/activate e si lancia

====================================================
from __future__ import absolute_import, division, print_function, unicode_literals

import numpy as np

import tensorflow as tf
assert tf.__version__.startswith('2')

from tensorflow_examples.lite.model_customization.core.data_util.image_dataloader import ImageClassifierDataLoader
from tensorflow_examples.lite.model_customization.core.task import image_classifier
from tensorflow_examples.lite.model_customization.core.task.model_spec import efficientnet_b0_spec
from tensorflow_examples.lite.model_customization.core.task.model_spec import ImageModelSpec

import matplotlib.pyplot as plt

image_path = "./foraminifera/"

data = ImageClassifierDataLoader.from_folder(image_path)
model = image_classifier.create(data)
loss, accuracy = model.evaluate()
model.export('fora_classifier.tflite', 'fora_labels.txt')
====================================================

a questo punto si e' pronti per testare la libreria con il seguente script
====================================================
"""label_image for tflite."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import argparse
import numpy as np

from PIL import Image

#import tensorflow as tf # TF2
import tflite_runtime.interpreter as tflite



def load_labels(filename):
  with open(filename, 'r') as f:
    return [line.strip() for line in f.readlines()]


if __name__ == '__main__':
  parser = argparse.ArgumentParser()
  parser.add_argument(
      '-i',
      '--image',
      default='./test_image/orbulina_test.jpg',
      help='image to be classified')
  parser.add_argument(
      '-m',
      '--model_file',
      default='fora_classifier.tflite',
      help='.tflite model to be executed')
  parser.add_argument(
      '-l',
      '--label_file',
      default='fora_labels.txt',
      help='name of file containing labels')
  parser.add_argument(
      '--input_mean',
      default=127.5, type=float,
      help='input_mean')
  parser.add_argument(
      '--input_std',
      default=127.5, type=float,
      help='input standard deviation')
  args = parser.parse_args()

  #interpreter = tf.lite.Interpreter(model_path=args.model_file)

  interpreter = tflite.Interpreter(model_path=args.model_file)

  interpreter.allocate_tensors()

  input_details = interpreter.get_input_details()
  output_details = interpreter.get_output_details()

  # check the type of the input tensor
  floating_model = input_details[0]['dtype'] == np.float32

  # NxHxWxC, H:1, W:2
  height = input_details[0]['shape'][1]
  width = input_details[0]['shape'][2]
  img = Image.open(args.image).resize((width, height))

  # add N dim
  input_data = np.expand_dims(img, axis=0)

  if floating_model:
    input_data = (np.float32(input_data) - args.input_mean) / args.input_std

  interpreter.set_tensor(input_details[0]['index'], input_data)

  interpreter.invoke()

  output_data = interpreter.get_tensor(output_details[0]['index'])
  results = np.squeeze(output_data)

  top_k = results.argsort()[-5:][::-1]
  labels = load_labels(args.label_file)
  for i in top_k:
    if floating_model:
      print('{:08.6f}: {}'.format(float(results[i]), labels[i]))
    else:
      print('{:08.6f}: {}'.format(float(results[i] / 255.0), labels[i]))
====================================================

lo script si lancia con la sintassi

(venv) luca@debian:~/tensor/tf_spectra/fora$ python classify_image_2.py -i ./test_image/crassaformis_test.jpg

in pratica con lo switch -i si indica il file immagine di test di cui si vuole riconoscere l'immagine

un risultato di esempio e'

(venv) luca@debian:~/tensor/tf_spectra/fora$ python classify_image_2.py -i ./test_image/crassaformis_test.jpg
INFO: Initialized TensorFlow Lite runtime.
0.823927: globorotalia_crassaformis
0.098882: globigerina_bulloides
0.031008: globigerinella_calida
0.023761: globiferina_falconensis
0.015816: orbulina_universa

in pratica l'immagine di test di una Globorotalia Crassaformis e' stato riconosciuta al 83% di confidenza con la giusta classificazione
Mettendo in tabella i risultati


In verde la corretta attribuzione, in arancione attribuzione errata da parte della rete neurale

commentando i risultati si osserva che su 12 test 8 hanno portato ad un corretto riconoscimento, con punteggi molto alti in caso di foraminiferi molto caratteristici come le Orbuline, in un caso e' stato individuato in modo corretto il genere ma non la specie, mentre nei due rimanenti casi di errata classificazione l'errore e' molto marcato

1. Hsiang AY, Brombacher A, Rillo MC, Mleneck-Vautravers MJ, Conn S, Lordsmith S, Jentzen A, Henehan MJ, Metcalfe B, Fenton I, Wade BS, Fox L, Meilland J, Davis CV, Baranowski U, Groeneveld J, Edgar KM, Movellan A, Aze T, Dowsett H, Miller G, Rios N, Hull PM. (2019) Endless Forams: >34,000 modern planktonic foraminiferal images for taxonomic training and automated species recognition using convolutional neural networks. Paleoceanography & Paleoclimatology, 34. https://doi.org/10.1029/2019PA003612
2. Elder L.E., Hsiang A.Y., Nelson K., Strotz L.C., Kahanamoku S.S., Hull P.M. Sixty-one thousand recent planktonic foraminifera from the Atlantic Ocean. Scientific Data 5: 180109. https://doi.org/10.1038/sdata.2018.109
3. Rillo M.C., Whittaker J., Ezard T.H.G., Purvis A., Henderson A.S., Stukins S., Miller C.G. 2016. The unknown planktonic foraminiferal pioneer Henry Buckley and his collection at The Natural History Museum, London. Journal of Micropalaeontology. https://doi.org/10.1144/jmpaleo2016-020

venerdì 4 maggio 2018

Intersezione superficie topografica con piani con qgSurf

Un plugin di QGis che non conoscevo ma che puo' tornare utile per verificare la correttezza delle faglie disegnate su una carta geologica.

Questo plugin, presi un modello digitale del terreno e un piano, calcola la linea di intersezione tra le due superfici...una regola della V tecnologica



Per tentativi puo' essere anche stimata la pendenza del piano di faglia ove sia riesca a mappare il contatto in campagna

sabato 7 aprile 2018

Test di Project Tango Tablet per la geologia

Scansioni 3D effettuate con Tablet Project Tango ed RTAB-MAP con visualizzazione in MeshLab

L'aspetto interessante di utilizzare RTAB-MAP e' che non e' necessario utilizzare una ripresa statica da cavelletto ma si puo' muovere la camera perche' viene gestito il motion tracking. Questo permette di riprendere anche zone che possono risultare "in ombra" da una ripresa frontale

Prova 1 







Prova 2







domenica 4 marzo 2018

Supervised Classification USGS Spectral Library

Dopo aver letto le note introduttive di TensorFlow (l'insieme di librerie di Machine Learning di Google) ed in particolar modo l'esempio sulla classificazione degli Iris sulla base di elementi morfologici mi e' venuta la domanda se il metodo era applicabile anche al telerilevamento iperspettrale, oggetto della mia tesi di dottorato.

Il problema a questo punto era prima di tutto trovare una base dati di training il piu' possibile popolata di elementi (in dottorato ho fatto collezione di spettri ma non organizzate in modo tale da essere utili allo scopo) e la scelta e' caduta su USGS Spectral Libray (gia' usata qui) in particolare per la sezione degli spettri del satellite Hyperion, oramai dismesso ma che e' stato anche lui oggetto di parte della tesi

Visto che a lavoro mi sto occupando in questi di minerali della famiglia dell'asbesto (http://debiaonoldcomputers.blogspot.it/2018/02/actinolite-e-tremolite-ad-impruneta.html) ho provato ad estrarre gli spettri di actinolite, tremolite e serpentino.
In totale il dataset, costituito dalle 3 classi, e' da 23 campioni da 234 bande (nel database gli spettri sono in numero maggiore ma sono stati scelti solo quelli che avevano un campionamento omogeneo)

Gli spettri sono stati tutti normalizzati prima di entrare nel file del dataset.
Si poteva a questo punto si fare l'analisi con le sole componenti principali (PCA analysis) per rendere il dataset piu' piccolo eliminando tutti i dati autocorrelati oppure dare in pasto all'algoritmo tutto lo spettro; ho provato con la seconda strada per vedere come si comportava il codice di calcolo

Nel dettaglio sono stati scelti

Actinolite
s07HYPRN_Actinolite_HS22.1B_ASDFRb_AREF
s07HYPRN_Actinolite_HS22.2B_ASDFRb_AREF
s07HYPRN_Actinolite_HS22.3B_ASDFRb_AREF
s07HYPRN_Actinolite_HS22.4B_ASDFRb_AREF
s07HYPRN_Actinolite_HS116.1B_ASDFRb_AREF
s07HYPRN_Actinolite_HS116.2B_ASDFRb_AREF
s07HYPRN_Actinolite_HS116.3B_ASDFRb_AREF
s07HYPRN_Actinolite_HS116.4B_ASDFRb_AREF
s07HYPRN_Actinolite_HS315.1B_ASDFRb_AREF
s07HYPRN_Actinolite_HS315.2B_ASDFRb_AREF

Tremolite

s07HYPRN_Tremolite_HS18.1B_ASDFRc_AREF
s07HYPRN_Tremolite_HS18.2B_ASDFRc_AREF
s07HYPRN_Tremolite_HS18.3_BECKc_AREF
s07HYPRN_Tremolite_HS18.3B_ASDFRc_AREF
s07HYPRN_Tremolite_HS18.3B_NIC4ccc_RREF
s07HYPRN_Tremolite_HS18.4B_ASDFRc_AREF
s07HYPRN_Tremolite_NMNH117611.HCl_BECKb_AREF
s07HYPRN_Tremolite_NMNH117611.HCL_NIC4bb_RREF

Serpentino
s07HYPRN_Serpentine_HS8.2B_ASDFRc_AREF
s07HYPRN_Serpentine_HS8.3B_ASDFRc_AREF
s07HYPRN_Serpentine_HS8.3B_BECKc_AREF
s07HYPRN_Serpentine_HS8.4B_ASDFRc_AREF
s07HYPRN_Serpentine_HS8.6_ASDFRc_AREF
s07HYPRN_Serpentine_HS318.1B_ASDFRc_AREF
s07HYPRN_Serpentine_HS318.2B_ASDFRc_AREF
s07HYPRN_Serpentine_HS318.3B_ASDFRc_AREF
s07HYPRN_Serpentine_HS318.4B_ASDFRc_AREF
s07HYPRN_Serpentine_HS318.4B_BECKc_AREF
s07HYPRN_Serpentine_HS318.6_ASDFRc_AREF



Ho provato ad usare TensorFlow per il calcolo ma mi e' decisamente ostico. Ho trovato che l'analisi del dataset Iris era stato trattato in questo post trattato con la libreria SkLearn,  con associate le librerie esterne Pandas e NumPy. I data sono tutti contenuti in un file, sia il training che il test dataset; il dataset sara' diviso in due dalla funzione train_test_split (in questo caso su 23 campioni viene scelto il 20% dei dati come parte di test del modello)

La prima riga del file all_data.txt e' costituito da un header con i nome delle 234 bande e un campo finale del nome (quindi la tabella e' di 235x24 celle)

-----------------------------------------
import pandas as pd
import numpy as np

dataset = pd.read_csv("all_data.txt")

X = dataset.iloc[:,:233].values
y = dataset['classe'].values

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.20, random_state = 82)

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)

# Fitting Multiclass Logistic Classification to the Training set
from sklearn.linear_model import LogisticRegression

logisticregression = LogisticRegression()
logisticregression.fit(X_train, y_train)

LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
verbose=0, warm_start=False)

# Predicting the Test set results
y_pred = logisticregression.predict(X_test)
print(y_pred)

#lets see the actual and predicted value side by side
y_compare = np.vstack((y_test,y_pred)).T

#actual value on the left side and predicted value on the right hand side
#printing the top 5 values
y_compare[:5,:]

# Making the Confusion Matrix
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(y_test, y_pred)
print(cm)

#finding accuracy from the confusion matrix.
a = cm.shape
corrPred = 0
falsePred = 0

for row in range(a[0]):
for c in range(a[1]):
if row == c:
corrPred +=cm[row,c]
else:
falsePred += cm[row,c]
print('Correct predictions: ', corrPred)
print('False predictions', falsePred)
print ('Accuracy of the multiclass logistic classification is: ', corrPred/(cm.sum()))
------------------------------

Il risultato finale e' seguente

[[2 0 0]
 [0 2 0]
 [0 0 1]]
('Correct predictions: ', 5)
('False predictions', 0)

('Accuracy of the multiclass logistic classification is: ', 1)


Per prova ho deliberatamente messo nel test dataset uno spettro con una classificazione sbagliata. Il risultato e' che l'algoritmo mi ha riportato un errore di classificazione...quindi direi che funziona

Il dataset e' sicuramente modesto e si possono fare miglioramenti ma direi che lo scopo e' stato raggiunto. La cosa interessante sarebbe ripetere questa esperienza con TensorFlow









mercoledì 21 febbraio 2018

Actinolite e Tremolite ad Impruneta

Avvertimento : l'actinolite appartiene alla famiglia dell'asbesto e quindi e' potenzialmente sorgente di fibre che possono generare malattie cancerogene quali il mesotelioma. Da gestire in condizioni di sicurezza 

Avvertimento : i dati sotto riportati sono relativi solo ad osservazioni. Non ho possibilita' di effettuare determinazione analitiche o al SEM per avere conferma della determinazione.
-----------------------------------------------------



Affioramento Sassi Neri Impruneta lungo la strada che collega a Strada in Chianti

L'actinolite e' un anfibolo  parente stretto della tremolite (cui si differenzia solo per la differente percentuale di sostituzione tra Mg e Fe). La formula completa, sia di tremolite che di actinolite, e'
Ca2(Mg,Fe2+)5Si8O22(OH)2

Gli anfiboli contenenti Ca sono tipici delle rocce ultrafemiche

Per i riferimenti di minerologia si rimanda a Mindat

A livello di affioramento l'unico indicatore e' dato dal colore che risulta essere verde tendente al nero per l'actinolite mentre e' bianco per la tremolite. Mentre nel campione 1 sembra che vi sia solo actinolite nel campione 2 potrebbe esservi sia tremolite che actinolite

Colpendo con il martello (cosa da non fare ma mi e' scappato un colpo) il colore diventa bianco

Tutte le mineralizzazioni si presentano con accrescimenti lungo la superficie di frattura e non ortogonali alla stessa.
Nel campione 1 l'aspetto dell'actinolite e' piu' massivo. Si nota la direzione di accrescimento ma non si notano in modo evidente la separazione tra i diversi cristalli lungo l'asse maggiore.

Campione 1

Campione 1

Campione 2
Alcuni ingrandimenti con un microscopio USB a 50X (circa). In affioramento con la lentina da 20x la struttura e' comunque evidente



Dettagli Campione 1







Dettagli Campione 2








Da notare che in alcune zone del campione al microscopio si nota una estesa presenza di fori nella zona esposta all'alterazione meteorica dell'affioramento...non ho idea se si tratti di azione di piante od animali

Sono presenti anche campioni di serpentino ma non ho visto presenza di crisotilo alla scala di affioramento

Campione 3






mercoledì 15 marzo 2017

mercoledì 8 febbraio 2017

Geological compass for Micro:bit

Aggiornamento qui
Un semplice programma in MicroPython per usare l'accelerometro e la bussola di Micro:bit come bussola da rilevamento geologico


Premendo il pulsante A si acquisisce la misura (se la bussola non e' stata calibrata sara' necessario farlo ruotando il dispositivo fino a creare un cerchio sul display). Il pulsante B ripete la visualizzazione dell'ultima misura
Il costo finale, compreso di batterie, e' di circa 25 euro



Il calcolo matematico di pitch e roll partendo dai dati dell'accelerometro e' stato ripreso da qui

Mark Pedley, Tilt Sensing Using a Three-AxisAccelerometer. Freescale Semiconductor Document Number: AN3461 Application Note Rev. 6, 03/2013


Il calcolo di strike e dip partendo da pitch e roll e' stato ripreso da questo articolo

 R.N. Barbosa *,1 , J.B. Wilkerson2 , H.P. Denton3 and D.C. Yoder 2, Slope Gradient and Vehicle Attitude Definition Based on Pitch and Roll Angle Measurements: A Simplified ApproachThe Open Agriculture Journal, 2012, 6, 36-40

in questa immagine l'orientazione degli assi dell'accelerometro e della bussola rispetto alla scheda (da notare che l'accelerazione Gz e' negativa perche' all'accelerometro e' saldato sulla parte opposta della scheda)




(nota : i valori di accelerazione vengono mostrati come numeri interi e non come float (m/s2) o come frazioni dell'accelerazione di gravita' (come in Android)
----------------------
from microbit import display
import microbit
import math

display.scroll("GEOCOMPASS")

cc = 180/math.pi
teta = 0
dip = 0

while True:
    if microbit.button_a.is_pressed():
        # if necessary calibrate the compass
        if not microbit.compass.is_calibrated():
            display.scroll("Calibrate compass")
            microbit.compass.calibrate()
            microbit.sleep(2500)
        gx = microbit.accelerometer.get_x()
        gy = microbit.accelerometer.get_y()
        gz = microbit.accelerometer.get_z()
        hd = microbit.compass.heading()
        # change the sign to uniform to Android
        gx = - gx
        gz = - gz
        pitch = math.atan2(gy, math.sqrt((gx*gx)+(gz*gz)))
        roll = math.atan2(-gx, gz)
        # from pitch/roll to strike/dip
        p2 = math.sin(pitch)*math.sin(pitch)
        r2 = math.sin(roll)*math.sin(roll)
        t1 = math.sqrt(p2+r2)
        # ---
        teta = (cc * math.asin(t1))
        sigma = math.asin(math.sin(roll)/t1)
        sigma = (cc * sigma)
        # -----
        # primo quadrante
        if ((gy <= 0) and (gx < 0)):
            sigma = sigma
        # secondo quadrante
        if ((gy > 0) and (gx < 0)):
            sigma = 180 - sigma
        # terzo quadrante
        if ((gy > 0) and (gx >= 0)):
            sigma = 180 - sigma
        # quarto quadrante
        if ((gy <= 0) and (gx >= 0)):
            sigma = 360 + sigma
        dip = (sigma + hd) % 360
        display.scroll("Dip " + str(int(teta)) + "/Str " + str(int(dip)))
        microbit.sleep(500)
    if microbit.button_b.is_pressed():
        # the B button repeat the last measure
        display.scroll("Dip " + str(int(teta)) + "/Str " + str(int(dip)))
----------------------

mercoledì 26 ottobre 2016

Frane a Pelago

Durante i periodo di tesi (e per le escursioni) mi sono trovato spesso a percorrere la strada che da Pontassieve porta alla Consuma e mi sono sempre chiesto il motivo della forma del versante indicato dalla freccia.


Da un estratto del Piano Strutturale del Comune di Pelago Tav. G.02 Novembre 2013 aggiornato a marzo 2014 l'area di interesse (al centro dell'immagine sottostante) e' circondato da orli di scarpata e corone di frana attive ed inattive


Guardando nella parte basale del pendio si vede quelli che sembrano gli accumuli al piede della frana

A questo punto mi sono incuriosito se si osservano movimenti dall'analisi dei permanent scatterers. I dati sono pubblici ma non facilmente reperibili. Si deve andare sul portale cartografico nazionale, selezionare il servizio WMS (il tema non e' listato in basso) e cercare Progetto Permanent Scatterer (http://wms.pcn.minambiente.it/pst) per poi selezionare la propria area (Firenze) ed il satellite di riferimento (ERS o CSK Cosmo SkyMed)


Da questa immagine non sembra che vi sia attivita' di versante



Da Peccioli ad Houston

In in precedente post avevo parlato della roccia lunare in esposizione a Peccioli mettendo l'accento sull'eta' dichiarata di 4 milioni di anni specialmente in relazione alla geologia lunare


del resto pero' avevo letto la scheda petrografica del campione (qui) che riportava la seguente tabella


e per quanto strano avevo accettato la datazione (potevano essere stati compiuti due errori dalla Nasa su due documenti di cui uno ufficiale??)

A distanza di 5 mesi ho ripreso in mano l'argomento e curiosando su internet ho trovato il sito geologialunare.it curato dal Dr. G. Turdo con cui ho avviato una corrispondenza via mail sull'argomento a cui ha gentilmente risposto (grazie mille) facendomi notare che nello stesso documento alla stessa pagina poco piu' in alto della tabella c'era il seguente grafico


in cui viene chiaramente indicata la data in b.y. (miliardi di anni) e non m.y (milioni di anni)...scemo io a non vederlo.

A questo punto ho scritto una mail al Lunar Sample Facility Laboratory che, in stile pienamente USA, mi ha risposto nel giro di un paio di giorni confermandomi l'errore e che la targhetta sarebbe stata corretta una volta che il campione rientrava nel loro laboratorio dopo le esposizioni a giro per  il mondo...concludendo con un invito a visitare il loro laboratorio ad Houston come ringraziamento

venerdì 21 ottobre 2016

Faglia bordiera del bacino di Firenze

Fin dai tempi dell'Universita' mi e' stata chiara la genesi del bacino di Firenze ma ultimamente mi sono reso conto che non sapevo ubicare effettivamente la posizione e non avevo nemmeno l'idea se ci fossero delle evidenze geomorfologiche....ho scoperto che ho sempre avuto sott'occhio questa struttura tettonica senza di fatto riconoscerla.

Per trovare la posizione della struttura, vista la sua dimensione, si parte semplicemente dai dati SRTM


Dopo averli banalmente elaborati in QGis si vede chiaramente la rottura di pendio che e' stata contraddistanta dalla linea blu. Sovrapponendo la foto area di Google Maps per avere dei riferimenti la situazione e' piu' o meno questa


Dato che abito in zona Fiesole ho zoomato nell'area per vedere se potevo riconoscere qualcosa di familiare. La rottura di pendio e' posta in prossimita' di S.Domenico


Andando sulla strada che da Firenze porta a Fiesole, fermandosi tra l'ospedale di Camerata e S.Domenico e guardando in direzione Est, si ossserva chiaramente l'andamento della struttura (da mangiarsi i gomiti....ci saro' passato mille volte)


Per conferma ho cercato qualcosa a riguardo in bibliografia. Qui e' evidenziato il sistema tettonico completo e si vede che la faglia risalta ancora di piu' perche' si trova al contatto tra i depositi Plio-Pleistocenici e il Macigno

Coli Rubellini 2007

Coli Rubellini 2013



Tour de France Firenze Gran Depart

  Poagcar Wout Van Aert Vingegaard       Van Der Poel