Projects/pr-iworms

From CIMENT
Jump to: navigation, search

Contents

[edit] pr-iworms: Innovative workflow for time synchronization of massive seismologic data

This part is a template getting data from Perseus. If you need to make some changes, please, go into the Perseus page of your project (link bellow).

[edit] Title/Titre: "Innovative workflow for time synchronization of massive seismologic data"

Perseus link / Lien Perseus: https://perseus.ujf-grenoble.fr/project/pr-iworms

[edit] Scientific description / Description scientifique

La sismologie moderne repose sur le déploiement de grands réseaux d’observation composés de milliers de capteurs indépendants enregistrant les données en continu. La qualité des données dépend de la synchronisation en quasi-temps réel des différentes stations qui composent ces réseaux. Celle-ci est cruciale pour l’étude des séismes et les problématiques d’imagerie et de monitoring. Les objectifs de ce projet interdisciplinaire sont (1) de mettre en œuvre une nouvelle méthode de correction de la datation basée sur l’inter-corrélation de bruit sismique ambiant (2) d’appliquer ce traitement à des flux continus de données.

[edit] Technical description / Description technique

preprocess python obspy séquentiel avec RAM négligeable corr + inv mtx Fortran openmp avec xx GB RAM (à voir selon taille du réseau)

luke 1 noeud pour les tests de prototypage io read sur montage nfs /media/resif io write sur scratch local de luke4 ou luke24/25

  • - Funding type: cnrs mastodons
  • - Parallelism: parallel
  • - Number of cores: max 1 node
  • - Memory: max 100GB/process
  • - Duration: 36
  • - Software needs: python avec obspy et h5py

openmp fftw lapack

  • - Visualization needs: non
  • - Storage needs:
  • - Training needs: non
  • - HPC-PME: NO
  • - CIGRI: YES
  • - Other needs:

[edit] Bibliography / Bibliographie

[edit] Scientific bibliography / Blibiliographie scientifique (EDIT)

@article{Sens_Sch_nfelder_2008,
  title = {Synchronizing seismic networks with ambient noise},
  volume = {174},
  issn = {1365-246X},
  url = {http://dx.doi.org/10.1111/j.1365-246X.2008.03842.x},
  doi = {10.1111/j.1365-246x.2008.03842.x},
  number = {3},
  journal = {Geophysical Journal International},
  publisher = {Oxford University Press (OUP)},
  author = {Sens-Schönfelder, Christoph},
  year = {2008},
  month = {sep},
  pages = {966–970}
}

[edit] Bibliography containing results obtained with Ciment / Bibliographie contenant des résultats obtenus grace à Ciment (EDIT)

No data. Please, fill from Perseus.

[edit] Activity report validation

When you have finished your activity report on the second part of this page (below), you have to send a notification to the manager of your CIMENT pole. To do so, you have to click on the provided button on the project page:

https://perseus.ujf-grenoble.fr/project/pr-iworms#activity-report

[edit] Bibliographie projet

[edit] 2017 activity report / Rapport d'activite 2017

[edit] Status

  • Date: 10/01/2018
  • Completion: 25%
  • Deadline extension needed: 2019-12-31

[edit] Results obtained

File:Mastodons2017 Trame Rapport Renouvellement Vfinale.pdf

[edit] 2018 activity report / Rapport d'activite 2018

[edit] Status

  • Date: 14/01/2019
  • Completion: 50%
  • Deadline extension needed: 2019-12-31

[edit] Results obtained

voir ci dessous les CR de réunions et sessions de travail ayant eu lieu en 2018

[edit] 2019 activity report / Rapport d'activite 2019

[edit] Status

  • Date: 13/01/2020
  • Completion: 90%
  • Deadline extension needed: 2020-12-31

[edit] Results obtained

[edit] Cahier de projet

[edit] Choix des données

Code source pour:

  • compter le nb de traces d'une liste donnée pour un jour donné dans /media/resif
  • rapatrier (rsync) cette liste de resif9@BIO vers luke:/nfs_scratch
$ svn co https://forge.osug.fr/svn/isterre-correlation/IWORMS/IWORMS_rapatrie_src IWORMS_rapatrie_src
A    IWORMS_rapatrie_src/rapatrie.sh
A    IWORMS_rapatrie_src/lookfornbstaJ.sh
A    IWORMS_rapatrie_src/stationsByCodes.txt
Checked out revision 201.

Réseau:

  • AlpArray: 121 stations -> on en choisit 50 sur critère d'interdistance minimum : zone de densité max: Alpes maritimes
  • 3 composantes
  • ~100Hz

Resif data portal

-> Seismic networks -> Browse by stations

-> Networks -> VIRTUAL -> ALPARRAY_FR

Search

Draw selection area

Search

Map -> Table

Export as text

#Network|Station|Location|Channel|Latitude|Longitude|Elevation|Depth|Azimuth|Dip|SensorDescription|Scale|ScaleFreq|ScaleUnits|SampleRate|StartTime|EndTime
FR|GRN|00|HHE|45.241300|5.743900|1040.0|0.0|90.0|0.0|Trillium Compact seismometer|2.9964e+08|1|M/S - meter per second (velocity)|100|2013-06-28T12:00:00.000000|2014-07-09T10:00:00.000000
FR|GRN|00|HHE|45.241300|5.743900|1040.0|0.0|90.0|0.0|Trillium Compact seismometer|2.9964e+08|1|M/S - meter per second (velocity)|100|2014-07-09T10:05:00.000000|2500-12-31T23:59:59.999999
FR|GRN|00|HHN|45.241300|5.743900|1040.0|0.0|0.0|0.0|Trillium Compact seismometer|2.9964e+08|1|M/S - meter per second (velocity)|100|2013-06-28T12:00:00.000000|2014-07-09T10:00:00.000000
FR|GRN|00|HHN|45.241300|5.743900|1040.0|0.0|0.0|0.0|Trillium Compact seismometer|2.9964e+08|1|M/S - meter per second (velocity)|100|2014-07-09T10:05:00.000000|2500-12-31T23:59:59.999999
...

plusieurs lignes pour une même trace car interruption dans la campagne d'acquisition (Attention ce n'est pas la même chose que des gaps !)

On récupère un fichier texte: stationsByCodes.txt

$ md5sum stationsByCodes.txt 
0a4a694904402509142da0821efef49a  stationsByCodes.txt
$ ll stationsByCodes.txt
-r--r--r-- 1 lecointre l-isterre 372491 May 24 14:40 stationsByCodes.txt

On définit une boite large avec 64 stations, donc 192 traces si on ne sélectionne que les HH[ENZ].

Box 64 stations AlpArray
(thumbnail)
Box 64 stations AlpArray
(thumbnail)
Box 64 stations AlpArray (zoom)


lecoinal@ist-oar:~/iWORMS_rapatrie_src$ sed 's/^M/\n/g' stationsByCodes.txt > zzzzz
lecoinal@ist-oar:~/iWORMS_rapatrie_src$ awk -F'|' '{print $1 " " $2 " " $3 " " $4}' zzzzz | grep -v "#Network" | grep "HH[ENZ]$"|sort -u|wc -l
216
lecoinal@ist-oar:~/IWORMS_rapatrie_src$ awk -F'|' '{print $1 " " $2 " " $3 " " $4}' zzzzz | grep -v "#Network" | grep "HH[ENZ]$"|awk '{print $1 " "$2}'|sort -u|wc -l
64
lecoinal@ist-oar:~/IWORMS_rapatrie_src$ awk -F'|' '{print $1 " " $2 " " $3 " " $4}' zzzzz | grep -v "#Network" | grep "HH[ENZ]$"|grep RUSF|sort -u
FR RUSF 00 HHE
FR RUSF 00 HHN
FR RUSF 00 HHZ
FR RUSF 01 HHE
FR RUSF 01 HHN
FR RUSF 01 HHZ
FR RUSF 03 HHE
FR RUSF 03 HHN
FR RUSF 03 HHZ
FR RUSF 04 HHE
FR RUSF 04 HHN
FR RUSF 04 HHZ
FR RUSF 05 HHE
FR RUSF 05 HHN
FR RUSF 05 HHZ
FR RUSF 06 HHE
FR RUSF 06 HHN
FR RUSF 06 HHZ
FR RUSF 07 HHE
FR RUSF 07 HHN
FR RUSF 07 HHZ

3*64 = 192 < 216 car plusieurs locations code pour station RUSF:

RUSF is moving.png
FR_RUSF HHE 43.929300 5.481900
FR_RUSF HHE 43.935700 5.483100
FR_RUSF HHE 43.936300 5.465600
FR_RUSF HHE 43.940300 5.474600
FR_RUSF HHE 43.941100 5.483900
FR_RUSF HHE 43.941200 5.483700
FR_RUSF HHN 43.929300 5.481900
FR_RUSF HHN 43.935700 5.483100
FR_RUSF HHN 43.936300 5.465600
FR_RUSF HHN 43.940300 5.474600
FR_RUSF HHN 43.941100 5.483900
FR_RUSF HHN 43.941200 5.483700
FR_RUSF HHZ 43.929300 5.481900
FR_RUSF HHZ 43.935700 5.483100
FR_RUSF HHZ 43.936300 5.465600
FR_RUSF HHZ 43.940300 5.474600
FR_RUSF HHZ 43.941100 5.483900
FR_RUSF HHZ 43.941200 5.483700
$ grep RUSF stationsByCodes_nl.txt | grep "HH[ENZ]"|awk -F'|' '{print $3 " " $4 " " $5 " " $6 " " $16 " " $17}'|sort -k 4 | grep 2500
03 HHE 43.936300 5.465600 2015-01-27T12:58:00.000000 2500-12-31T23:59:59.999999
03 HHN 43.936300 5.465600 2015-01-27T12:58:00.000000 2500-12-31T23:59:59.999999
03 HHZ 43.936300 5.465600 2015-01-27T12:58:00.000000 2500-12-31T23:59:59.999999
04 HHE 43.940300 5.474600 2015-10-13T08:13:00.000000 2500-12-31T23:59:59.999999
04 HHN 43.940300 5.474600 2015-10-13T08:13:00.000000 2500-12-31T23:59:59.999999
04 HHZ 43.940300 5.474600 2015-10-13T08:13:00.000000 2500-12-31T23:59:59.999999
06 HHE 43.929300 5.481900 2015-01-27T10:37:00.000000 2500-12-31T23:59:59.999999
06 HHN 43.929300 5.481900 2015-01-27T10:37:00.000000 2500-12-31T23:59:59.999999
06 HHZ 43.929300 5.481900 2015-01-27T10:37:00.000000 2500-12-31T23:59:59.999999
05 HHE 43.935700 5.483100 2015-01-27T11:32:00.000000 2500-12-31T23:59:59.999999
05 HHN 43.935700 5.483100 2015-01-27T11:32:00.000000 2500-12-31T23:59:59.999999
05 HHZ 43.935700 5.483100 2015-01-27T11:32:00.000000 2500-12-31T23:59:59.999999
01 HHE 43.941200 5.483700 2015-01-27T15:20:00.000000 2500-12-31T23:59:59.999999
01 HHN 43.941200 5.483700 2015-01-27T15:20:00.000000 2500-12-31T23:59:59.999999
01 HHZ 43.941200 5.483700 2015-01-27T15:20:00.000000 2500-12-31T23:59:59.999999
07 HHE 43.941100 5.483900 2015-01-27T13:03:00.000000 2500-12-31T23:59:59.999999
07 HHN 43.941100 5.483900 2015-01-27T13:03:00.000000 2500-12-31T23:59:59.999999
07 HHZ 43.941100 5.483900 2015-01-27T13:03:00.000000 2500-12-31T23:59:59.999999
$

On considérera qu'un trace est identifiée par RESEAU.STATION.LOCATION.CHANNEL et qu'une trace ne peut pas changer de position (latitude,longitude) au cours du temps.

On cherche quelques jours consécutifs avec au moins 150 traces valides sur les 216:

lecoinal@ist-oar:~/IWORMS_rapatrie_src$ for j in $(seq -w 1 200) ; do echo -n "nb traces in /media/resif/validated_seismic_data : 2017$j : " ; ./lookfornbstaJ.sh 2017$j ; done
nb traces in /media/resif/validated_seismic_data : 2017001 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017002 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017003 : 168 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017004 : 168 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017005 : 168 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017006 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017007 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017008 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017009 : 171 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017010 : 168 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017011 : 171 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017012 : 171 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017013 : 171 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017014 : 171 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017015 : 171 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017016 : 168 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017017 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017018 : 165 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017019 : 168 zfilesfrom.txt
nb traces in /media/resif/validated_seismic_data : 2017020 : 168 zfilesfrom.txt
...

On choisit les 5 jours 2017011 - 2017015 (171 traces).

On rapatrie ces données sur luke:/nfs_scratch

lecoinal@ist-oar:~/iWORMS_rapatrie_src$ for jd in $(seq -w 10 16) ; do ./rapatrie.sh 20170$jd ; done > zzzzz

Box 72 stations avec composantes HH[ENZ]. On en a 57 qui fonctionnent sur la période 2017011-2017015 (donc 171 séries temporelles):

Box 72 stations AlpArray
(thumbnail)
Box 72 stations AlpArray

EDIT juillet 2017: crash de la baie MD1200 luke4 (donc perte des données prétraitées) et down du serveur luke4 (donc plus d'accès nfs au /media/resif resif9@BIO).

Il faut recommencer le preprocess sur un autre noeud luke : note luke24 accède désormais aussi à /media/resif.

Note: 3 nouvelles traces sont apparues au catalogue...

lecoinal@ist-oar:~/IWORMS_rapatrie_src$ date
Tue Jul 18 10:51:12 CEST 2017
lecoinal@ist-oar:~/IWORMS_rapatrie_src$ for j in $(seq -w 10 10) ; do echo -n "nb traces in /media/resif/validated_seismic_data : 20170$j : " ; ./lookfornbstaJ.sh 20170$j ; done
nb traces in /media/resif/validated_seismic_data : 2017010 : 171 zfilesfrom.txt
2.35018157958984375000 GB
                   encoding: STEIM 1 Compression (val:10)
nb traces in /media/resif/validated_seismic_data : 2017011 : 174 zfilesfrom.txt
2.35340118408203125000 GB
                   encoding: STEIM 1 Compression (val:10)
nb traces in /media/resif/validated_seismic_data : 2017012 : 174 zfilesfrom.txt
2.34509277343750000000 GB
                   encoding: STEIM 1 Compression (val:10)
nb traces in /media/resif/validated_seismic_data : 2017013 : 174 zfilesfrom.txt
2.50910568237304687500 GB
                   encoding: STEIM 1 Compression (val:10)
nb traces in /media/resif/validated_seismic_data : 2017014 : 174 zfilesfrom.txt
2.47692108154296875000 GB
                   encoding: STEIM 1 Compression (val:10)
nb traces in /media/resif/validated_seismic_data : 2017015 : 174 zfilesfrom.txt
2.40940093994140625000 GB
                   encoding: STEIM 1 Compression (val:10)
nb traces in /media/resif/validated_seismic_data : 2017016 : 171 zfilesfrom.txt
2.42947769165039062500 GB
                   encoding: STEIM 1 Compression (val:10)

[edit] Lien resif (RESIF9@BIO) - CIMENT

Demande d'export /media/resif en NFS sur luke4, luke24, luke25: Pierre Volcke, Bruno Bzeznik, Romain Cavagna : en cours...

En attendant: rsync via script rapatrie.sh (voir ci dessus)

  • AL: je souhaiterais exporter /media/resif vers luke4, luke24 et luke25.
  • PV: peux-tu confirmer qu'il s'agit bien de : luke4.u-ga.fr luke24.u-ga.fr luke25.u-ga.fr ? Sur l'export NFS en question, les données sont lisibles par le groupe POSIX dont le GID est 2136 (ceci est contraint par le fait que ces données doivent être lues par un groupe d'usagers sur le cluster ISterre, qui détient ce GID). Y a t-il moyen, du côté CIMENT, de faire en sorte que ce soit lisible aussi ? (peut être en faisant un remappage d'UID ?).
  • BB: Luke[4,24,25] sont bien des noeuds de l'équipe Ondes. Le remappage vers le gid d'un projet ciment (10343(pr-iworms)) serait certainement le plus propre. Je ne suis pas expert en remapping, mais ce remappage n'est-il pas à faire côté serveur NFS avec une entrée export du genre:
/media/resif luke4.ujf-grenoble.fr(rw,all_squash,anongid=10343)
  • PV: j'ai fait l'export:
[root@resif9 ~]# exportfs -v
...
/mnt/local/vg01/lv00
luke4.u-ga.fr(ro,async,wdelay,root_squash,all_squash,no_subtree_check,anonuid=99,anongid=10343,sec=sys,ro,root_squash,all_squash)
...

l'export s'appelle localement /mnt/local/vg01/lv00 mais on le nomme "traditionnellement" /media/resif/ du côté client. peux-tu tester sur luke4, et si ça marche je propage pour le reste ?

  • BB: Le montage se fait et l'usager "lecointre" peut copier des fichiers!
lecointre@luke4:~$ cp
/resif/validated_seismic_data/2015/7H/PBOU/DHZ.D/7H.PBOU.06.DHZ.D.2015.174 /tmp
lecointre@luke4:~$

Le problème c'est que j'ai l'impression que tout le monde en fait à accès aux fichiers...

bzizou@luke4:~$ cp
/resif/validated_seismic_data/2015/7H/PBOU/DHZ.D/7H.PBOU.06.DHZ.D.2015.174 /tmp
bzizou@luke4:~$

Et on voit:

bzizou@luke4:~$ ls -lhd /resif/validated_seismic_data/2015
drwxr-s--- 13 99 2136 104 Jul  7  2016 /resif/validated_seismic_data/2015

Donc, on dirait que le mapping n'a aucun effet, mais que les droits sur les répertoire non plus! Le "all_squash" est peut être en trop?

  • PV: merci pour le test, j'ai modifié les parmètres en retirant le all_squash :
luke4.uga.fr
(ro,async,wdelay,root_squash,no_subtree_check,anonuid=99,anongid=10343,sec=sys,ro,root_squash,no_all_squash)

ça dit quoi ?

  • BB: Zut, maintenant, même Albanne ne voit plus les fichiers :-(
bzizou@luke4:~$ ls -lhd /resif/validated_seismic_data/2015
ls: cannot access /resif/validated_seismic_data/2015: Permission denied
bzizou@luke4:~$ logout
root@luke4:~# su - lecointre
lecointre@luke4:~$ ls -lhd /resif/validated_seismic_data/2015
ls: cannot access /resif/validated_seismic_data/2015: Permission denied
lecointre@luke4:~$
root@luke4:~# ls -lhd /resif/validated_seismic_data/2015
drwxr-s--- 13 99 2136 104 Jul  7  2016 /resif/validated_seismic_data/2015

Pas évident ce mapping d'uid visiblement. J'ai vu sur le net passer des trucs avec le fichier /etc/idmapd.conf côté client... ça me semble bizarre. Bon, pour ne pas perdre de temps, solution crado, mais efficace:

# Creation du groupe resif en local
root@luke4:~# groupadd -g 2136 resif
# Ajout de l'utilisateur lecointre dans le groupe resif
root@luke4:~# vigr
# Test
lecointre@luke4:~$ ls -lhd /resif/validated_seismic_data/2015
drwxr-s--- 13 99 resif 104 Jul  7  2016 /resif/validated_seismic_data/2015

Albanne, peux-tu nous dire si cette solution convient pour toi?

  • AL : depuis luke4 je vois bien /resif (pas depuis luke24 ni luke25 par contre, mais dans l'immédiat je propose de laisser tomber ces 2 noeuds supplémentaires, les premiers proto pour iworms laissent penser que les calculs vont être assez légers, un noeud suffirait largement). Je peux donc maintenant laisser tomber de mon workflow cette étape (rapatrie rsync d'une liste de fichiers depuis ist-oar:/resif vers luke4:/nfs_scratch)  :
$ svn co https://forge.osug.fr/svn/isterre-correlation/IWORMS/IWORMS_rapatrie_src IWORMS_rapatrie_src
A    IWORMS_rapatrie_src/rapatrie.sh

et passer direct à l'étape suivante du workflow (preprocess) :

$ svn co https://forge.osug.fr/svn/isterre-correlation/IWORMS/IWORMS_preprocess_src IWORMS_preprocess_src
A    IWORMS_preprocess_src/submit_IWORMS_preprocess.sh
A    IWORMS_preprocess_src/whitenWhisper.py
A    IWORMS_preprocess_src/stationsByCodes.txt
A    IWORMS_preprocess_src/iWORMS_preprocess.py
Révision 202 extraite.

et changer juste une ligne du submit

- file="/nfs_scratch/lecointre/IWORMS/media/resif/validated_seismic_data/$YEAR/$net/$sta/$chan.D/$net.$sta.$loc.$chan.D.$YEAR.$JDAY"
+ file="/resif/validated_seismic_data/$YEAR/$net/$sta/$chan.D/$net.$sta.$loc.$chan.D.$YEAR.$JDAY"

et je vérifie que je regénère bien un fichier preprocessé identique au précédent :

# le précédent (input rapatriee au préalable sur /nfs_scratch)
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$ ll 2017.011.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 5615091 Jun 13 13:13 2017.011.trace.int8.h5
# le nouveau
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$ ll /nfs_scratch/lecointre/IWORMS_preprocess_3308363/2017.011.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 5615091 Jun 14 2017 /nfs_scratch/lecointre/IWORMS_preprocess_3308363/2017.011.trace.int8.h5
# ils contiennent la même liste de dataset :
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$ h5ls 2017.011.trace.int8.h5 > list1
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$ h5ls /nfs_scratch/lecointre/IWORMS_preprocess_3308363/2017.011.trace.int8.h5 > list2
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$ diff list1 list2
# et les datasets ont les mêmes valeurs
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$ h5diff /nfs_scratch/lecointre/IWORMS_preprocess_3308363/2017.011.trace.int8.h5 2017.011.trace.int8.h5
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2$

A noter que les perfs de preprocessing ne semblent pas impactées par le fait de lire les inputs dans luke:/nfs_scratch ou dans resif9.ujf-grenoble.fr:/mnt/local/vg01/lv00. Donc ca convient.

  • PV : Je me demande si, dans la conf serveur, la dernière ligne (anongid=10343) est vraiment prise en compte ou si la première (anongid=2136) prend le pas (malgré ce que prétend exportfs -v).
/mnt/local/vg01/lv00 \
         -ro,async,no_subtree_check,anonuid=99,anongid=2136 \
        ist-calcul21.univ-savoie.fr \
         ist-calcul22.univ-savoie.fr \
        ist-calcul23.univ-savoie.fr
         luke4.ujf-grenoble.fr(ro,anongid=10343) <===

Pour les autres montages sur les calculateurs ISterre, Rodolphe et Jean-Noel n'ont pas eu, à ma connaissance, eu besoin de jouer avec la conf de leur client NFS. En conclusion, si cela te convient je repasse à :

/mnt/local/vg01/lv00 \
         -ro,async,no_subtree_check,anonuid=99,anongid=2136 \
        ist-calcul21.univ-savoie.fr \
         ist-calcul22.univ-savoie.fr \
        ist-calcul23.univ-savoie.fr
         luke4.ujf-grenoble.fr <===

...compte tenu du fait que tu as crée le groupe 2136 localement. Par ailleurs il se peut que d'autres usagers CIMENT soit un jour intéressé pour accèder à /media/resif : tu souhaiteras alors peut être régler plus proprement ce problème de conf.

  • BB : Jean-Noel, tu as un mapping de gid qui fonctionne, ou est-ce que tu as créé un groupe local avec un gid identique comme j'ai fini par faire?
  • JNB : pour ce qui est du mapping entre CIMENT et ISTerre voici ce que nous faisons : nous creons dans notre LDAP un groupe cime-pro-{projet} ou cime-equ-{equipe} dont le gid est le meme que le gid CIMENT. Rien de special dans les options de montage NFS. Dois-je comprendre que le projet pr-iworms est destine a ce type d'usage ?
  • BB : OK, je vois, donc en fait, tu as résolu le problème inverse. C'est à dire que tu calque ton gid local sur le gid du projet Ciment. Là, c'est resif qui nous impose un gid... Merci pour l'info!

Oui, sauf que le gid côté résif est déja fixé à 2136 et qu'il ne correspond pas au gid du projet i-worms (10343)

  • PV : pas tout à fait : le GID (2136) du côté du serveur NFS portant les données RESIF correspond au GID déclaré sur les serveurs de calcul ISterre pour le groupe ayant accès à ces données (groupe resif-data), (...) c'est pour cela que nous avons tenté, dans notre conf, de faire le remapping d'UID pour coller à l'UID 10343 (CIMENT), sans pour autant impacter l'export vers les serveurs ISterre. Or, il semble que ce

remapping ne fonctionne pas (pour situer Jean-Noel : c'est là où nous en étions au moment où tu es rentré dans la boucle).

@Bruno : est ce que le "hack" que tu as mis en place (déclaration en local du groupe 2136 sur CIMENT) te satisfait ? (=> est-ce acceptable d'en rester là dans la mesure où Albanne peut désormais accèder aux datas ?) Si non : sur la machine RESIF nous n'avons aucune contrainte particulière sur l'UID et le GID. Y a t-il un moyen d'utiliser le même GID pour le groupe ISterre et le groupe CIMENT ? D'autres idées ?

  • BB : Disons que c'est moyennement satisfaisant. L'idéal serait que quelqu'un devienne expert des problématiques de mapping de gid en NFS, j'espère qu'il y a des solutions plus élégantes ;-)

Maintenant, ça fonctionne et ne pose pas de problème de maintenance particulier dans l'immédiat, donc cela ne me gène pas.

  • JNB : @Pierre : tu peux aussi inverser le point de vue et :
- creer un groupe avec id 10343
- chgroup -R 10343 /resif
- me le faire savoir pour que je cree le groupe 10343 dans LDAP OSUG
- me dire qui je dois mettre dans le groupe

A toi de voir.

  • PV : ok. Je propose un second essai avec ces mapping : J'ai ré-écrit la table d'exports de façon un peu différente (j'avais remarqué par le passé que l'ordre de déclaration pouvait compter... sans

trouver d'explication.).

@Bruno : si cette nouvelle configuration fonctionne, le groupe 2136 que tu as déclaré localement ne sera plus nécessaire. Peux-tu le retirer, et vérifier à nouveau si le compte d'Albanne a accès depuis son groupe CIMENT originel, et que les autres users n'ont pas accès ?

Si cet essai n'a pas plus de succès que le précèdent, une question sur la suggestion de Jean-Noël :

oui je peux passer tout ce répertoire sous le GID 10343,
mais si je comprends bien ce sont donc tous les users du groupe 2136 (resif-dat) que je te demanderais alors de basculer dans le groupe 10343 (afin que l'accès aux datas soit préservés depuis les machines ISterre) ?
  • JNB : oui
  • BB : Je pense que ça ne marche pas mieux, je vois toujours les fichiers dans le groupe 2136 :
lecointre@luke4:~$ ls -ln /resif/validated_seismic_data/2017/Z3/ |tail
drwxr-xr-x 5 99 2136 42 Feb  1 19:26 A208A
drwxr-xr-x 5 99 2136 42 May 30 10:56 A209B
drwxr-xr-x 5 99 2136 42 Jan  9 19:43 A210A
drwxr-xr-x 5 99 2136 42 Jan  6 18:05 A211A
drwxr-xr-x 5 99 2136 42 Jan  3 18:03 A212A
drwxr-xr-x 5 99 2136 42 Jan  9 19:43 A213A
drwxr-xr-x 5 99 2136 42 Jan  3 18:03 A214A
drwxr-xr-x 5 99 2136 42 Jan  3 18:03 A215A
drwxr-xr-x 5 99 2136 42 Jan  5 18:03 A216A
drwxr-xr-x 5 99 2136 42 Jan  3 18:03 A217A
  • PV : Oui, les fichiers appartiennent effectivement à 99.2136. Ce qu'il faut vérifier, c'est si Albanne a bien à nouveau accès en lecture sans pour autant être dans le groupe 2136.

L'inconvénient est (comme mentionné par Bruno dans un précèdent mail) que tous les users ayant accès à ce montage NFS sur la machine Luke auront aussi accès aux datas... En fait je ne pense pas que l'on puisse faire plus fin avec ce simple remapping NFS, mais je ne suis pas expert. Pour revenir sur la proposition de Jean-Noel (basculer toutes ces données sous le groupe 10343 en créant un groupe LDAP adéquat de notre côté) : on aurait donc un GID homogène entre CIMENT et OSUG/ISterre, problème réglé. Mais : comme CIMENT gère des GID par projet (10343(pr-iworms)), je suppose que l'on va retrouver la même problématique dès qu'un autre projet CIMENT souhaitera accéder aux mêmes datas ?

  • BB : Je viens de faire un test en enlevant Albanne du groupe 2136 et elle accède toujours aux fichiers. Mais effectivement, le problème, c'est que tout le monde accède aux fichiers... j'ai fait un test avec mon login (bzizou).
root@luke4:~# id lecointre
uid=10057(lecointre) gid=10015(l-isterre)
groups=10015(l-isterre),10028(m-froggy),10045(li-maimosine),10044(li-users),10047(m-r2d2),10051(m-gofree),10053(m-wiki),10059(m-cigri),10064(pr-formation-hpc-130614),10101(pr-beamforming),10057(li-cigri),10117(li-luke),10126(m-luke),10146(pr-whisper),10172(pr-sanjacinto),10174(pr-imvort),10186(pr-formation-ced-mdm-2015),10201(pr-vi-hps-tw18),10236(pr-eventdetection),10245(pr-rsc),10267(p-page-astro-geo),10305(pr-imagin),10307(f-ondes),10318(pr-formation-ced-2017),10343(pr-iworms),10348(pr-metaforet)
root@luke4:~# id lecointre |grep resif
root@luke4:~# su - lecointre
lecointre@luke4:~$ ls -l
/resif/validated_seismic_data/2017/FR/ARBF/HHZ.D/FR.ARBF.00.HHZ.D.2017.011
-r--r--r-- 1 99 2136 23318528 Jan 18 04:30
/resif/validated_seismic_data/2017/FR/ARBF/HHZ.D/FR.ARBF.00.HHZ.D.2017.011
lecointre@luke4:~$ logout
root@luke4:~# su - bzizou
bzizou@luke4:~$ ls -l
/resif/validated_seismic_data/2017/FR/ARBF/HHZ.D/FR.ARBF.00.HHZ.D.2017.011
-r--r--r-- 1 99 2136 23318528 Jan 18 04:30
/resif/validated_seismic_data/2017/FR/ARBF/HHZ.D/FR.ARBF.00.HHZ.D.2017.011
bzizou@luke4:~$

Je remet albanne dans le groupe 2136 en attendant...

  • PV : un bilan des manips des derniers jours, pour mémoire il y avait plusieurs points à traiter:

- faire en sorte que les données soient accessibles à la fois depuis les machines ISterre et depuis les machines CIMENT, en tenant compte de la différence des référentiels utilisateur.

- tenir compte du fait que les usagers peuvent appartenir à plus d'un groupe projet (au sens POSIX), notamment dans CIMENT (=> les droits POSIX classiques ne sont pas suffisants pour exprimer cela).

- conserver la "privacy" des données.

Jusqu'à présent, l'approche avait été d'utiliser la fonction de remapping d'UID/GID dans NFS, mais cette technique n'a pas fonctionné de façon satisfaisante du côté CIMENT. Ce n'était de toute façon qu'une rustine, dont on aurait vite atteint les limites.

Bruno a suggéré de mettre en place des droits complémentaires sous formes d'ACL => c'est ce que nous avons fait, et supprimé le remapping d'UID. Cela a fonctionné côté CIMENT, mais par effet de bord nous avons perdu les accès côté ISterre. Suite à échange avec Rodolphe, le problème a été localisé dans la conf du serveur NFS sur resif9 (=> correction faite).

Donc : les modalités d'accès restent les mêmes depuis le cluster ISterre (groupe LDAP dédié) et les accès aux usagers CIMENT sont gérés par ajout d'ACL, au cas par cas. Cela reste donc souple. En pj une doc pour la gestion de cette configuration, si vous voulez bien jeter un coup d'oeil (https://wiki-geodata.obs.ujf-grenoble.fr/doku.php?id=geodata:systemes_et_reseaux:inventaire_serveurs:resif9:media_resif).

Bruno : peux-tu nous rappeller le GID du groupe dédié à iWorms ? Albanne : à ma connaissance les données sont visibles depuis luke4 et luke24, cela peut être étendu à d'autres noeuds.

[edit] InternalMSEEDReadingWarning lors de la lecture (python obspy) de données requêtées sur RESIF@SUMMER

[Voir ticket ouvert chez RESIF-dc]

lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$ grep InternalMSEEDReadingWarning OAR.cigri.10970.*.stderr
OAR.cigri.10970.3745480.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=121711, xn=-11308
OAR.cigri.10970.3745480.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 10739712 is not valid SEED. The rest of the file will not be read.
OAR.cigri.10970.3759430.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=122199, xn=-10572
OAR.cigri.10970.3759430.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 10727424 is not valid SEED. The rest of the file will not be read.
OAR.cigri.10970.3759431.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=63421, xn=-9986
OAR.cigri.10970.3759431.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 10883072 is not valid SEED. The rest of the file will not be read.
OAR.cigri.10970.3759432.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=125123, xn=-7986
OAR.cigri.10970.3759432.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 10383360 is not valid SEED. The rest of the file will not be read.
OAR.cigri.10970.3759577.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=93949, xn=96685
OAR.cigri.10970.3761064.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=17043, xn=-10642
OAR.cigri.10970.3761064.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 10592256 is not valid SEED. The rest of the file will not be read.
OAR.cigri.10970.3761065.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=182305, xn=-10560
OAR.cigri.10970.3761065.stderr:/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 8781824 is not valid SEED. The rest of the file will not be read.
lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$

Il s'agit de jobs ayant tourné sur la période 14-17 novembre 2017:

lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$ for j in $(grep InternalMSEEDReadingWarning OAR.cigri.10970.*.stderr|awk -F. '{print $4}'|sort -u) ; do echo "--- $j ---" ; oarstat -fj $j | grep Time ; oarstat -fj $j | grep hostname ; done
--- 3745480 ---
   submissionTime = 2017-11-14 19:53:14
   startTime = 2017-11-14 19:53:19
   stopTime = 2017-11-14 20:00:56
   assigned_hostnames = luke8
--- 3759430 ---
   submissionTime = 2017-11-16 13:28:20
   startTime = 2017-11-16 13:28:31
   stopTime = 2017-11-16 13:41:10
   assigned_hostnames = luke30
--- 3759431 ---
   submissionTime = 2017-11-16 13:28:20
   startTime = 2017-11-16 13:28:31
   stopTime = 2017-11-16 13:40:37
   assigned_hostnames = luke29
--- 3759432 ---
   submissionTime = 2017-11-16 13:28:20
   startTime = 2017-11-16 13:28:31
   stopTime = 2017-11-16 13:40:54
   assigned_hostnames = luke29
--- 3759577 ---
   submissionTime = 2017-11-16 14:28:13
   startTime = 2017-11-16 14:28:25
   stopTime = 2017-11-16 14:41:38
   assigned_hostnames = luke31
--- 3761064 ---
   submissionTime = 2017-11-17 02:20:43
   startTime = 2017-11-17 02:20:45
   stopTime = 2017-11-17 02:26:45
   assigned_hostnames = luke8
--- 3761065 ---
   submissionTime = 2017-11-17 02:20:43
   startTime = 2017-11-17 02:20:45
   stopTime = 2017-11-17 02:26:45
   assigned_hostnames = luke8
lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$

On a traité 366 jours (1 an de l'été 2016 à l'été 2017) et à peu près 78 traces par jour, on a traité un total de 25,874 traces exactement, et ces Warning ne se sont produits que pour 7 traces journalières.

lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$ for f in $(ls OAR.cigri.10970.*.stdout) ; do tag=$(head -1 $f) ; echo -n "$tag " ; head -5 $f | tail -1 | awk  '{print $8}'; done|sort -u|awk '{print $2}'|sed "/^$/d"|wc -l
366
lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$ for f in $(ls OAR.cigri.10970.*.stdout) ; do tag=$(head -1 $f) ; echo -n "$tag " ; head -5 $f | tail -1 | awk  '{print $8}'; done|sort -u|awk '{print $2}'|sed "/^$/d"|awk "{SUM+=$1}END{print SUM}"
25874

Ces 7 traces sont les suivantes:

1) 6 traces pour lesquelles on a du 2 type de Warning:

  • InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=<x>, xn=<y>
  • InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset <z> is not valid SEED. The rest of the file will not be read.

2) 1 trace pour laquelle on a eu un unique type de Warning:

  • InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=<x>, xn=<y>

Voir sujets connexes (?) sur github obspy: https://github.com/obspy/obspy/issues/962 https://github.com/obspy/obspy/issues/1981


1) La liste des 6 traces journalières en cause est:

  • 3745480 network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000 :
  • 3759430 network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-11T23:59:00.000000&endtime=2017-03-13T00:01:00.000000
  • 3759431 network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-12T23:59:00.000000&endtime=2017-03-14T00:01:00.000000
  • 3759432 network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-13T23:59:00.000000&endtime=2017-03-15T00:01:00.000000
  • 3761064 network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-05-16T23:59:00.000000&endtime=2017-05-18T00:01:00.000000
  • 3761065 network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-05-17T23:59:00.000000&endtime=2017-05-19T00:01:00.000000

Exemple détaillé d'une trace pour laquelle on a les deux types de warning: FR.ARBF.00.HHZ 2016-10-04T23:59:00.000000 - 2016-10-06T00:01:00.000000:

converted 'http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000' (ANSI_X3.4-1968) -> 'http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000' (UTF-8)
--2017-11-14 19:51:05--  http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 401
Reusing existing connection to ws.resif.fr:80.
HTTP request sent, awaiting response... 200
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

    0K .......... .......... .......... .......... .......... 91.8M
   50K .......... .......... .......... .......... .......... 2.42K
...
 10350K .......... .......... .......... .......... ..........  209M
10400K .......... .......... .......... .......... ..........  184M
10450K .......... .......... .......... ........               191M=21s

2017-11-14 19:51:26 (505 KB/s) - Read error at byte 10739759 (Success).Retrying.

--2017-11-14 19:51:27--  (try: 2)  http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 401
Reusing existing connection to ws.resif.fr:80.
HTTP request sent, awaiting response... 200
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

    0K .......... .......... .......... .......... .......... 65.7M
   50K .......... .......... .......... .......... .......... 61.9M
  100K .......... .......... .......... .......... .......... 70.1M
...
22750K .......... .......... .......... .......... ..........  163M
22800K .......... .......... .......... .......... ..........  137M
22850K .......... .......... .......... .......... ..........  178M
22900K .......... ..                                           129M=0.2s

2017-11-14 19:51:27 (114 MB/s) - written to stdout [23461888]

/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=121711, xn=-11308
 warnings.warn(*_i)
/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Record starting at offset 10739712 is not valid SEED. The rest of the file will not be read.
 warnings.warn(*_i)

On constate qu'il y a eu deux essais successifs de requête, avec à chaque fois une erreur autour de la position 10739712 ou 10739759.

Les positions correspondantes pour les 6 traces détaillées ci dessus sont:

lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$ grep "Read error at byte" *err
OAR.cigri.10970.3745480.stderr:2017-11-14 19:51:26 (505 KB/s) - Read error at byte 10739759 (Success).Retrying.
OAR.cigri.10970.3759430.stderr:2017-11-16 13:27:57 (421 KB/s) - Read error at byte 10727471 (Success).Retrying.
OAR.cigri.10970.3759431.stderr:2017-11-16 13:28:02 (427 KB/s) - Read error at byte 10883127 (Success).Retrying.
OAR.cigri.10970.3759432.stderr:2017-11-16 13:28:02 (409 KB/s) - Read error at byte 10383407 (Success).Retrying.
OAR.cigri.10970.3761064.stderr:2017-11-17 02:18:43 (489 KB/s) - Read error at byte 10592319 (Success).Retrying.
OAR.cigri.10970.3761065.stderr:2017-11-17 02:18:43 (405 KB/s) - Read error at byte 8781863 (Success).Retrying.

2) La trace pour laquelle on a eu un unique warning est:

  • network=FR&station=OGDI&location=00&channel=HHZ&starttime=2017-03-16T23:59:00.000000&endtime=2017-03-18T00:01:00.000000
converted 'http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=OGDI&location=00&channel=HHZ&starttime=2017-03-16T23:59:00.000000&endtime=2017-03-18T00:01:00.000000' (ANSI_X3.4-1968) -> 'http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=OGDI&location=00&channel=HHZ&starttime=2017-03-16T23:59:00.000000&endtime=2017-03-18T00:01:00.000000' (UTF-8)
--2017-11-16 14:33:59--  http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=OGDI&location=00&channel=HHZ&starttime=2017-03-16T23:59:00.000000&endtime=2017-03-18T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 401
Reusing existing connection to ws.resif.fr:80.
HTTP request sent, awaiting response... 200
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

    0K .......... .......... .......... .......... ..........  332K
   50K .......... .......... .......... .......... .......... 26.7K
  100K .......... .......... .......... .......... ..........  139M
...
18200K .......... .......... .......... .......... .......... 42.6M
18250K .......... .......... .......... .......... .......... 38.2M
18300K .......... .......... .......... ..                    40.3M=2.2s

2017-11-16 14:34:07 (8.17 MB/s) - written to stdout [18771968]

/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=93949, xn=96685
 warnings.warn(*_i)

Donc on n'a pas eu les deux essais successifs de requête.

1) Concernant les 6 traces avec double Warning:

  • Aujourd'hui je ne parviens pas à reproduire ces warning (observés en nov 2017).
  • Les outputs 1bit après prétraitement de ces données aujourd'hui diffère de celui effectué en nov 2017 sur ces "mêmes" (?) données (environ 3 à 4 % de +/-1 ont changé de signe)
    • Est ce que ca signifie qu'il y a eu un update de ces données depuis nov 2017 ?

Exemple: au 21/06/2018 les warning ne se reproduisent plus:

$ source /applis/ciment/v2/env.bash 
$ module load python/2.7.12_gcc-4.6.2
$ cat readobspy.py 
#!/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.6/bin/python
# -*- coding: utf-8 -*-
## cat file.mseed.gz | gunzip | python2.7 readobspy.py
import sys
import obspy.core
from StringIO import StringIO
stringio_obj = StringIO(sys.stdin.read())
st=obspy.core.read(stringio_obj)
print st
$ wget 'http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000' -O - | python2.7 readobspy.py 
--2018-06-21 13:53:06--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2016-10-04T23:59:00.000000&endtime=2016-10-06T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [  <=>                                                                                                    ]  22.38M  4.21MB/s    in 5.3s    

2018-06-21 13:53:11 (4.21 MB/s) - written to stdout [23461888]

3 Trace(s) in Stream:
FR.ARBF.00.HHZ | 2016-10-04T23:58:50.119400Z - 2016-10-05T00:00:00.607400Z | 125.0 Hz, 8812 samples
FR.ARBF.00.HHZ | 2016-10-05T00:00:00.664700Z - 2016-10-06T00:00:00.656700Z | 125.0 Hz, 10800000 samples
FR.ARBF.00.HHZ | 2016-10-06T00:00:00.714200Z - 2016-10-06T00:01:01.122200Z | 125.0 Hz, 7552 samples
$ url='http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-11T23:59:00.000000&endtime=2017-03-13T00:01:00.000000'
$ wget $url -O - | python2.7 readobspy.py
--2018-06-21 13:56:49--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-11T23:59:00.000000&endtime=2017-03-13T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [     <=>                                                                                                 ]  22.28M  4.37MB/s    in 5.1s    

2018-06-21 13:56:59 (4.37 MB/s) - written to stdout [23363584]

3 Trace(s) in Stream:
FR.ARBF.00.HHZ | 2017-03-11T23:58:56.872000Z - 2017-03-12T00:00:00.480000Z | 125.0 Hz, 7952 samples
FR.ARBF.00.HHZ | 2017-03-12T00:00:00.538800Z - 2017-03-13T00:00:00.530800Z | 125.0 Hz, 10800000 samples
FR.ARBF.00.HHZ | 2017-03-13T00:00:00.589500Z - 2017-03-13T00:01:01.445500Z | 125.0 Hz, 7608 samples
$ url='http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-12T23:59:00.000000&endtime=2017-03-14T00:01:00.000000'
$ wget $url -O - | python2.7 readobspy.py 
--2018-06-21 13:57:35--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-12T23:59:00.000000&endtime=2017-03-14T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [     <=>                                                                                                 ]  22.29M  4.12MB/s    in 5.4s    

2018-06-21 13:57:45 (4.12 MB/s) - written to stdout [23371776]

3 Trace(s) in Stream:
FR.ARBF.00.HHZ | 2017-03-12T23:58:44.954800Z - 2017-03-13T00:00:00.530800Z | 125.0 Hz, 9448 samples
FR.ARBF.00.HHZ | 2017-03-13T00:00:00.589500Z - 2017-03-14T00:00:00.581500Z | 125.0 Hz, 10800000 samples
FR.ARBF.00.HHZ | 2017-03-14T00:00:00.640100Z - 2017-03-14T00:01:01.976100Z | 125.0 Hz, 7668 samples
$ url='http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-13T23:59:00.000000&endtime=2017-03-15T00:01:00.000000'
$ wget $url -O - | python2.7 readobspy.py
--2018-06-21 13:58:09--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-03-13T23:59:00.000000&endtime=2017-03-15T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [     <=>                                                                                                 ]  22.29M  4.17MB/s    in 5.3s    

2018-06-21 13:58:19 (4.17 MB/s) - written to stdout [23367680]

3 Trace(s) in Stream:
FR.ARBF.00.HHZ | 2017-03-13T23:58:51.949500Z - 2017-03-14T00:00:00.581500Z | 125.0 Hz, 8580 samples
FR.ARBF.00.HHZ | 2017-03-14T00:00:00.640100Z - 2017-03-15T00:00:00.632100Z | 125.0 Hz, 10800000 samples
FR.ARBF.00.HHZ | 2017-03-15T00:00:00.690800Z - 2017-03-15T00:01:01.658800Z | 125.0 Hz, 7622 samples
$ url='http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-05-16T23:59:00.000000&endtime=2017-05-18T00:01:00.000000'
$ wget $url -O - | python2.7 readobspy.py
--2018-06-21 13:58:39--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-05-16T23:59:00.000000&endtime=2017-05-18T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [     <=>                                                                                                 ]  22.06M  4.18MB/s    in 5.3s    

2018-06-21 13:58:49 (4.18 MB/s) - written to stdout [23134208]

3 Trace(s) in Stream:
FR.ARBF.00.HHZ | 2017-05-16T23:58:47.401900Z - 2017-05-17T00:00:00.833900Z | 125.0 Hz, 9180 samples
FR.ARBF.00.HHZ | 2017-05-17T00:00:00.893200Z - 2017-05-18T00:00:00.885200Z | 125.0 Hz, 10800000 samples
FR.ARBF.00.HHZ | 2017-05-18T00:00:00.944400Z - 2017-05-18T00:01:02.568400Z | 125.0 Hz, 7704 samples
$ url='http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-05-17T23:59:00.000000&endtime=2017-05-19T00:01:00.000000'
$ wget $url -O - | python2.7 readobspy.py
--2018-06-21 13:59:04--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-05-17T23:59:00.000000&endtime=2017-05-19T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [    <=>                                                                                                  ]  22.07M  4.27MB/s    in 5.2s    

2018-06-21 13:59:14 (4.27 MB/s) - written to stdout [23146496]

3 Trace(s) in Stream:
FR.ARBF.00.HHZ | 2017-05-17T23:58:57.133200Z - 2017-05-18T00:00:00.885200Z | 125.0 Hz, 7970 samples
FR.ARBF.00.HHZ | 2017-05-18T00:00:00.944400Z - 2017-05-19T00:00:00.936400Z | 125.0 Hz, 10800000 samples
FR.ARBF.00.HHZ | 2017-05-19T00:00:00.995500Z - 2017-05-19T00:01:03.019500Z | 125.0 Hz, 7754 samples
$

Exemple: au 21/06/2018 on ne reproduit pas les outputs de nov 2017: 3% (=283,872/8,640,000) de données 1bit diffèrent.

$ h5diff 10970/IWORMS_preprocess_10970_3745480_2016279/2016.279.trace.int8.h5 11743/IWORMS_preprocess_11743_7573524_2016279/2016.279.trace.int8.h5 /FR.ARBF.00.HHZ
dataset: </FR.ARBF.00.HHZ> and </FR.ARBF.00.HHZ>
283872 differences found

2) Concernant la trace avec simple warning, j'arrive à le reproduire encore aujourd'hui:

$ url='http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=OGDI&location=00&channel=HHZ&starttime=2017-03-16T23:59:00.000000&endtime=2017-03-18T00:01:00.000000'
$ wget $url -O - | python2.7 readobspy.py
--2018-06-21 14:07:40--  http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=OGDI&location=00&channel=HHZ&starttime=2017-03-16T23:59:00.000000&endtime=2017-03-18T00:01:00.000000
Resolving ws.resif.fr (ws.resif.fr)... 152.77.1.122
Connecting to ws.resif.fr (ws.resif.fr)|152.77.1.122|:80... connected.
HTTP request sent, awaiting response... 200 
Length: unspecified [application/vnd.fdsn.mseed]
Saving to: 'STDOUT'

-                                                    [  <=>                                                                                                    ]  17.90M  3.31MB/s    in 5.4s    

2018-06-21 14:07:46 (3.31 MB/s) - written to stdout [18771968]

/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=93949, xn=96685
 warnings.warn(*_i)
2 Trace(s) in Stream:
FR.OGDI.00.HHZ | 2017-03-16T23:58:51.065064Z - 2017-03-16T23:59:58.995064Z | 100.0 Hz, 6794 samples
FR.OGDI.00.HHZ | 2017-03-17T00:02:10.625058Z - 2017-03-17T23:59:58.995058Z | 100.0 Hz, 8626838 samples
$ 

Et j'obtiens le même message d'erreur en remplacant la requête wget par python.obspy.fdsn.Client:

$ cat iWORMS_preprocess1.py
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
from obspy.core import read,UTCDateTime
import sys
cyr = sys.argv[1]  # as a string : "%04d"
cjd = sys.argv[2]  # as a string : "%03d"
yr = int(cyr)
jd = int(cjd)
opt = int(sys.argv[3])
if (opt == 0) or (opt == 3): 
	# Read stream from input files from NFS /media/resif (opt=0) or from input stream from wget -O - (opt=3)
	#s = read('tmp.mseed')
	from StringIO import StringIO
	try:
		s = read(StringIO(sys.stdin.read()))
	except Exception as e:
		print e
		sys.exit(0)  # cigri should not abort
elif opt == -1:
	# Used only for debug
	s = read("tmp.mseed")
elif opt == 1:
        # Download input stream thanks to obspy.clients.fdsn
	inet=sys.argv[4]
	ista=sys.argv[5]
	iloc=sys.argv[6]
	ichan=sys.argv[7]
	from obspy.clients.fdsn import Client
	import ConfigParser
	config=ConfigParser.ConfigParser()
	import getpass
	USERNAME=getpass.getuser()
	CONFIG_FILE="/home/"+USERNAME+"/.iworms.conf"
	import os
	os.path.isfile(CONFIG_FILE)
	config.read(CONFIG_FILE)
	USER=config.get('global','user')
	PASS=config.get('global','password')
	client = Client("RESIF",user=USER,password=PASS)
	t = UTCDateTime(cyr+cjd+"T00:00:00.000")
	try:
		s = client.get_waveforms(inet, ista, iloc, ichan, t - 60, t + 86400 + 60) # get from J-1min to J+1min to prevent bad time vector...
	except Exception as e:
		print e
		sys.exit(0)  # cigri should not abort
print(s)
$ python2.7 iWORMS_preprocess1.py 2017 076 1 FR OGDI 00 HHZ
/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=93949, xn=96685
 warnings.warn(*_i)
2 Trace(s) in Stream:
FR.OGDI.00.HHZ | 2017-03-16T23:58:51.065064Z - 2017-03-16T23:59:58.995064Z | 100.0 Hz, 6794 samples
FR.OGDI.00.HHZ | 2017-03-17T00:02:10.625058Z - 2017-03-17T23:59:58.995058Z | 100.0 Hz, 8626838 samples
$

Et la trace prétraitée (1bit) aujourd'hui (campagne 11743) est rigoureuseuement identique à celle prétraitée en nov 2017 (campagne 10970):

$ h5diff -v 10970/IWORMS_preprocess_10970_3759577_2017076/2017.076.trace.int8.h5 11743/IWORMS_preprocess_11743_7573528_2017076/2017.076.trace.int8.h5 /FR.OGDI.00.HHZ
dataset: </FR.OGDI.00.HHZ> and </FR.OGDI.00.HHZ>
0 differences found

[edit] Schéma workflow

Schema-2.png

[edit] Preprocessing des données

[edit] Methode 1

  • cleaning time vector:
    • concatenate J-1,J,J+1
    • for each trace: interpolate (linear) to 100Hz with an eventual new starttime : mod(starttime,100Hz)=0
    • cut between 00:00:00.000000 and 23:59:59.990000
  • highpass filter 0.05Hz (over the day if no gap or over each chunk of the day if gaps)
  • gapfill: fill gap with 0. TODO: keep the info about gaps position.
  • bandpass: 0.2-0.3333Hz (=3-5sec) (over the day)
  • decimation: 100Hz -> 10Hz (1sample/10 as for Japan, no spline interp)
  • whitening: 0.2-0.3333Hz (over the day)
  • check for NaN values
  • 1bit
  • Int8 encoding and rename output data channel with respect to channel naming SEED convention (HH[ENZ] -> BH[ENZ])
  •  write output as a dataset into HDF5 daily file with gzip compression

[edit] Methode 2

cf méthode L. Stehly (à vérifier avec lui)

  • cleaning time vector:
    • concatenate J-1,J,J+1
    • for each trace: interpolate (linear) to 100Hz with an eventual new starttime : mod(starttime,100Hz)=0
    • cut between 00:00:00.000000 and 23:59:59.990000
  • highpass filter 0.05Hz (over the day if no gap or over each chunk of the day if gaps)
  • gapfill: fill gap with 0. TODO: keep the info about gaps position.
  • SUM over several frequency bands:
    • bandpass: 0.3333-1Hz (1-3sec) + 0.2-0.3333Hz (=3-5sec) + 0.1429-0.2Hz (=5-7sec) + 0.1-0.1429Hz (=7-10sec)
    • normalisation par l'enveloppe (norm(trace/abs(hilbert(trace))))
  • decimation: 100Hz -> 10Hz (1sample/10 as for Japan, no spline interp)
  • no whitening
  • check for NaN values
  • 1bit
  • Int8 encoding and rename output data channel with respect to channel naming SEED convention (HH[ENZ] -> BH[ENZ])
  • write output as a dataset into HDF5 daily file with gzip compression

[edit] Brouillon

Pb1 : Pourquoi ARBF est en 125Hz ?

Pb2 : Il faut recaler sur des centièmes de secondes pleines ? Est ce qu'on sait si ce décalage est fixe ou pas ? apparemment oui mais uniquement pour un segment de trace donné. mais pas forcément pour deux segments consécutifs (voir AJAC entre 2015 et 2016...). Et ce décalage n'est pas forcément fixe pour les trois canaux d'une même station (voir ARTF)

Pb3 : ca déborde sur J+1 ? Est il possible que ca déborde aussi sur J-1 ?

Pb4 : Il faut parfois aller chercher le début de J dans J-1 ? on ne peut pas travailler indépendamment d'une journée à l'autre. il faudra télécharger 3 jours pour pretraiter un jour (IOread *3)

lecointre@luke:~/iWORMS_preprocess_src$ msi -tg /nfs_scratch/lecointre/IWORMS/2016/FR/*/*/*
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHE    2016,001,00:00:01.913162 2016,002,00:00:00.303162  ==  100 8639840
FR_AJAC_00_HHN    2016,001,00:00:01.393162 2016,002,00:00:01.213162  ==  100 8639983
FR_AJAC_00_HHZ    2016,001,00:00:00.193161 2016,002,00:00:01.063161  ==  100 8640088
FR_ARBF_00_HHE    2016,001,00:00:00.180500 2016,002,00:00:00.172500  ==  125 10800000
FR_ARBF_00_HHN    2016,001,00:00:00.180500 2016,002,00:00:00.172500  ==  125 10800000
FR_ARBF_00_HHZ    2016,001,00:00:00.180500 2016,002,00:00:00.172500  ==  125 10800000
FR_ARTF_00_HHE    2016,001,00:00:01.420660 2016,002,00:00:01.370660  ==  100 8639996
FR_ARTF_00_HHN    2016,001,00:00:01.040659 2016,002,00:00:00.430659  ==  100 8639940
FR_ARTF_00_HHZ    2016,001,00:00:02.030659 2016,002,00:00:00.800659  ==  100 8639878
FR_ASEAF_00_HHE   2016,001,00:00:01.300000 2016,002,00:00:00.670000  ==  100 8639938
FR_ASEAF_00_HHN   2016,001,00:00:01.940000 2016,002,00:00:01.570000  ==  100 8639964
FR_ASEAF_00_HHZ   2016,001,00:00:02.140000 2016,002,00:00:01.490000  ==  100 8639936
FR_ATE_00_HHE     2016,001,00:00:03.115534 2016,001,08:38:03.395534  ==  100 3108029
FR_ATE_00_HHE     2016,001,08:38:03.435534 2016,002,00:00:01.685534 0.04 100 5531826
FR_ATE_00_HHN     2016,001,00:00:01.545534 2016,001,08:38:03.045534  ==  100 3108151
FR_ATE_00_HHN     2016,001,08:38:03.085534 2016,002,00:00:04.035534 0.04 100 5532096
FR_ATE_00_HHZ     2016,001,00:00:02.585534 2016,001,08:38:04.765534  ==  100 3108219
FR_ATE_00_HHZ     2016,001,08:38:04.805534 2016,002,00:00:01.795534 0.04 100 5531700
Total: 15 trace(s) with 18 segment(s)
lecointre@luke:~/iWORMS_preprocess_src$
lecoinal@ist-oar:/media/resif/validated_seismic_data$ msi -tg 2016/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2016.001 
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHZ    2016,001,00:00:00.193161 2016,002,00:00:01.063161  ==  100 8640088
Total: 1 trace(s) with 1 segment(s)
lecoinal@ist-oar:/media/resif/validated_seismic_data$ msi -tg 2015/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2015.365 
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHZ    2015,365,00:00:01.093163 2016,001,00:00:00.183163  ==  100 8639910
Total: 1 trace(s) with 1 segment(s)
lecoinal@ist-oar:/media/resif/validated_seismic_data$

Solution 1 : dataselect -Sd -A

lecoinal@ist-oar:/media/resif/validated_seismic_data$ dataselect -Sd 2015/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2015.365 2016/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2016.001 -A ~/tmp/C/%Y/%n/%s/%c.%q/%n.%s.%l.%c.%q.%Y.%j
lecoinal@ist-oar:/media/resif/validated_seismic_data$ ll ~/tmp/C/201*/*/*/*/*
-rw-r--r-- 1 lecoinal iste-equ-ondes 19488768 May 15 11:38 /home/lecoinal/tmp/C/2015/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2015.365
-rw-r--r-- 1 lecoinal iste-equ-ondes 19132416 May 15 11:38 /home/lecoinal/tmp/C/2016/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2016.001
-rw-r--r-- 1 lecoinal iste-equ-ondes     4096 May 15 11:38 /home/lecoinal/tmp/C/2016/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2016.002
lecoinal@ist-oar:/media/resif/validated_seismic_data$ msi -tg /home/lecoinal/tmp/C/2016/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2016.001
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHZ    2016,001,00:00:00.003163 2016,001,23:59:59.993161  ==  100 8640000
Total: 1 trace(s) with 1 segment(s)
lecoinal@ist-oar:/media/resif/validated_seismic_data$ 

Mais il faudrait réinterpoler sur des centièmes de sec entières auparavant ?

Solution 2 : dataselect -ts -te -A

lecoinal@ist-oar:/media/resif/validated_seismic_data$ dataselect  -ts 2016,001,00,00,00,000000 -te 2016,001,23,59,59,990000 2015/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2015.365 2016/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2016.001 -A ~/tmp/D/%Y/%n/%s/%c.%q/%n.%s.%l.%c.%q.%Y.%j
lecoinal@ist-oar:/media/resif/validated_seismic_data$ ll ~/tmp/D/201*/*/*/*/*
-rw-r--r-- 1 lecoinal iste-equ-ondes     4096 May 15 11:41 /home/lecoinal/tmp/D/2015/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2015.365
-rw-r--r-- 1 lecoinal iste-equ-ondes 19128320 May 15 11:41 /home/lecoinal/tmp/D/2016/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2016.001
lecoinal@ist-oar:/media/resif/validated_seismic_data$ msi -tg /home/lecoinal/tmp/D/2016/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2016.001
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHZ    2016,001,00:00:00.193161 2016,002,00:00:01.063161  ==  100 8640088
Total: 1 trace(s) with 1 segment(s)
lecoinal@ist-oar:/media/resif/validated_seismic_data$ msi -tg /home/lecoinal/tmp/D/2015/FR/AJAC/HHZ.M/FR.AJAC.00.HHZ.M.2015.365
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHZ    2015,365,23:59:58.093163 2016,001,00:00:00.183163  ==  100 210
Total: 1 trace(s) with 1 segment(s)
lecoinal@ist-oar:/media/resif/validated_seismic_data$ 

C'est pas bon.

Solution 3 : dataselect -ts -te -o

lecoinal@ist-oar:/media/resif/validated_seismic_data$ dataselect  -ts 2016,001,00,00,00,000000 -te 2016,001,23,59,59,990000 2015/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2015.365 2016/FR/AJAC/HHZ.D/FR.AJAC.00.HHZ.D.2016.001 -o ~/tmp/tmp.mseed
lecoinal@ist-oar:/media/resif/validated_seismic_data$ msi -tg ~/tmp/tmp.mseed
   Source                Start sample             End sample        Gap  Hz  Samples
FR_AJAC_00_HHZ    2015,365,23:59:58.093163 2016,002,00:00:01.063161  ==  100 8640298
Total: 1 trace(s) with 1 segment(s)

Une solution ?:

  • 0. cat inputfileJm1 inputfileJ [inputfileJp1] > tmp.mseed
  • Est ce que cette étape est utile ? 1. dataselect -ts YYYY,JJJ,00,00,00,000000 -te YYYY,JJJ,23,59,59,990000 tmp.mseed -o tmp2.mseed -> on obtient un seul fichier avec J (et des petites bouts de J-1 et J+1)
  • 2. obspy -> resample sur des centièmes de sec précises: tr.interpolate(sampling_rate=100.,starttime=UTCDateTime(year=$YYYY, julday=$JJJ),method='weighted_average_slopes')
lecoinal@ist-oar:/media/resif/validated_seismic_data$ python2.7
Python 2.7.9 (default, Jun 29 2016, 13:08:31) 
[GCC 4.9.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from obspy.core import read,UTCDateTime
>>> st=read('/home/lecoinal/tmp/tmp2.mseed')
>>> st.interpolate(sampling_rate=100.,starttime=UTCDateTime(year=2016, julday=1),npts=8640000,method='weighted_average_slopes')
<obspy.core.stream.Stream object at 0x7f77a003cd10>
>>> print st
1 Trace(s) in Stream:
FR.AJAC.00.HHZ | 2016-01-01T00:00:00.000000Z - 2016-01-01T23:59:59.990000Z | 100.0 Hz, 8640000 samples
>>> 
  • 3. dataselect -Sd -A : inutile si option npts de obspy

[edit] Source code

Seule la méthode 1 est codée pour l'instant

$ svn co https://forge.osug.fr/svn/isterre-correlation/IWORMS/IWORMS_preprocess_src IWORMS_preprocess_src
A    IWORMS_preprocess_src/submit_IWORMS_preprocess.sh
A    IWORMS_preprocess_src/whitenWhisper.py
A    IWORMS_preprocess_src/stationsByCodes.txt
A    IWORMS_preprocess_src/iWORMS_preprocess.py
Révision 202 extraite.

[edit] Computing

Soumission preprocess jd 2017011 à 2017015:

$ oarsub -S "./submit_IWORMS_preprocess.sh 2017011"
$ oarsub -S "./submit_IWORMS_preprocess.sh 2017012"
$ oarsub -S "./submit_IWORMS_preprocess.sh 2017013"
$ oarsub -S "./submit_IWORMS_preprocess.sh 2017014"
$ oarsub -S "./submit_IWORMS_preprocess.sh 2017015"

ToDo:

  • output dataset name (STA_CHAN) : add network and location !!!
  • -> cigri
  • param: only J (J-1 and J+1 should be automatically updated)

Performance:

  • luke4
  • IO (read and write) on /nfs_scratch
  • ~5.3sec / 1day trace (if filestream input from NFS /resif)
  • 171 valid traces / day -> 0.25 CPUh (or elapsed h) / day
  • memory : (time -f "\t%M Maximum resident set size (KB)") : between 588MB and 1.9GB (depending on the trace)

Available options for getting input datasets:

  • stream input from NFS /resif (opt=0)
  • download input stream thanks to obspy.clients.fdsn (opt=1)
  • wget POST method (all traces in the same stream) (opt=2) : NON FONCTIONNEL
  • stream input from wget -O - (opt=3)

Options 1 and 3 give slightly different 1Bit results from option 0.

Get Input + Preprocess 2017011, 216 traces
PreprocessIWORMS NFS WS.png

Question: Ces 3 traces :

G SSB 10 HHE : 142 / 216 : no data .001604427 sec          |    G SSB 10 HHE : 142 / 216 : 3.246117723 sec
G SSB 10 HHN : 143 / 216 : no data .001365180 sec          |    G SSB 10 HHN : 143 / 216 : 3.228650985 sec
G SSB 10 HHZ : 144 / 216 : no data .001405344 sec          |    G SSB 10 HHZ : 144 / 216 : 3.267849534 sec

pour la journée 2017011, sont accessibles par webservices (wget method GET ou ospy.fdsn.client) mais ne sont pas présentes dans NFS /media/resif. pourquoi ?

Output:

$ ls -lh 2017.011.trace.int8.h5 
-rw-r--r-- 1 lecointre l-isterre 6.9M Jun  7 19:48 2017.011.trace.int8.h5
$ h5ls -v 2017.011.trace.int8.h5
Opened "2017.011.trace.int8.h5" with sec2 driver.
FR.ARBF.00.BHE           Dataset {864000/864000}
    Location:  1:800
    Links:     1
    Chunks:    {13500} 13500 bytes
    Storage:   864000 logical bytes, 39758 allocated bytes, 2173.15% utilization
    Filter-0:  fletcher32-3  {}
    Filter-1:  deflate-1 OPT {6}
    Type:      native signed char
FR.ARBF.00.BHN           Dataset {864000/864000}
    Location:  1:43254
    Links:     1
    Chunks:    {13500} 13500 bytes
    Storage:   864000 logical bytes, 39976 allocated bytes, 2161.30% utilization
    Filter-0:  fletcher32-3  {}
    Filter-1:  deflate-1 OPT {6}
    Type:      native signed char
...

Deux bandes de fréquence (bandpass filter):

  • 3-5sec (0.2-0.3333Hz):
luke4:~$ ls -h /scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.2_0.3333
-rw-r--r-- 1 lecointre l-isterre 6.9M Jun  7 19:48 2017.011.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 6.9M Jun  7 19:49 2017.012.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 6.9M Jun  7 19:49 2017.013.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 6.9M Jun  7 19:49 2017.014.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 7.0M Jun  7 19:50 2017.015.trace.int8.h5
  • 5-10sec (0.1-0.2Hz):
luke4:~$ ls -h /scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2
-rw-r--r-- 1 lecointre l-isterre 5.4M Jun 13 13:13 2017.011.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 5.4M Jun 13 13:13 2017.012.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 5.4M Jun 13 13:14 2017.013.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 5.4M Jun 13 13:13 2017.014.trace.int8.h5
-rw-r--r-- 1 lecointre l-isterre 5.4M Jun 13 13:14 2017.015.trace.int8.h5

[edit] Discussion: Bande de fréquence pour le preprocess et Coherence des fonctions de Green

[edit] Corrélations

Tests sur 2 bandes de fréquence:

  • 3-5sec : 0.2-0.3333Hz
  • 5-10sec : 0.1-0.2Hz
2017011 - calcul de correlations toutes les 30min (no overlap) - lag 300sec à 10Hz
EXP pair distance 3-5sec  : 5-10sec
EXP0 FR.GRN.00.BHZ FR.OGS2.00.BHZ 20182.705
EXP0 20km.png
EXP0 FR.GRN.00.BHZ FR.OGAG.00.BHZ 80430.430
EXP0 80km.png
EXP0 FR.OGGM.00.BHZ FR.RUSF.03.BHZ 150441.797
EXP0 150km.png
EXP0 FR.SPIF.00.BHZ Z3.A215A.00.BHZ 220310.312
EXP0 220km.png
EXP0 Z3.A176A.00.BHZ Z3.A204A.00.BHZ 301677.375
EXP0 300km.png
EXP1 FR.RUSF.01.BHZ FR.RUSF.03.BHZ 3114.258
EXP1 3km.png
EXP2 FR.SAOF.00.BHZ FR.TURF.00.BHZ 13027.878
EXP2 13km.png
EXP3 Z3.A205A.00.BHZ Z3.A206A.00.BHZ 29870.939
EXP3 30km.png

[edit] Cohérence des corrélations de GFs

Plot (pour ces bandes de fréquence), pour chaque paire, et en fonction de la distance (en km):

mean5d ( maxm ( ABS ( corrlag5sec (GFm,lag300sec , dailystackGFm,lag300sec) ) ) )

avec m=30min ; 60min ; 120min

1 seule itération

Pour le plot on garde la valeur absolue de cette cohérence avant de moyenner sur 5 jours.

2017011:2017015 - 5day mean absolute value of coherence between corr(30/60/120min) and daily corrstack - max corr in lags +/-5sec
3-5sec (0.2-0.3333Hz) 5-10sec (0.1-0.2Hz) 10-20sec (0.05-0.1Hz)
30min (f(distance))
BF COH30.3.5.png
BF COH30.5.10.png
BF COH30.10.20.png
30min (distribution)
BF COH30.3.5 hist.png
BF COH30.5.10 hist.png
BF COH30.10.20 hist.png
60min
BF COH60.3.5.png
BF COH60.5.10.png
BF COH60.10.20.png
60min (distribution)
BF COH60.3.5 hist.png
BF COH60.5.10 hist.png
BF COH60.10.20 hist.png
120min
BF COH120.3.5.png
BF COH120.5.10.png
BF COH120.10.20.png
120min (distribution)
BF COH120.3.5 hist.png
BF COH120.5.10 hist.png
BF COH120.10.20 hist.png


2017011:2017015 - All 5day coherence values between GF(30/60/120min) and daily GFstack - max abs corr in lags +/-5sec
Distribution des cohérences et approx loi normale Distribution des cohérences positives et négatives (Attention binning différents +/-) et approx loi normale Distribution de abs(cohérences) et approx loi normale 3-5sec (0.2-0.3333Hz) 5-10sec (0.1-0.2Hz) 10-20sec (0.05-0.1Hz)
30min
NCOH30. coh hist.png
COH30. posnegcoh hist.png
COH30. abscoh hist.png
NCOH30.3.5.png
NCOH30.5.10.png
NCOH30.10.20.png
60min
NCOH60. coh hist.png
COH60. posnegcoh hist.png
COH60. abscoh hist.png
NCOH60.3.5.png
NCOH60.5.10.png
NCOH60.10.20.png
120min
NCOH120 coh hist.png
COH120 posnegcoh hist.png
COH120 abscoh hist.png
NCOH120.3.5.png
NCOH120.5.10.png
NCOH120.10.20.png


[edit] Suivi d'une expérience: 5-10sec, 60min

  • 2017012
  • maxm ( ABS ( corrlag5sec (GFm,lag300sec , dailystackGFm,lag300sec) ) ), avec m60min
  • on garde ABS(COH) pour le plot
  • bande de frequence: 5-10sec
  • 1H-correlations
  • 2 pairs
    • coherence élevée: FR.OGGM.00.BHZ FR.OGSM.00.BHZ : 55.62 km : mean24h(coh_2017012) = 0.7563 (itération1) puis 0.7565 (itération2) 0.6496
    • coherence faible: FR.OGCB.00.BHN Z3.A192A.00.BHE : 55.62 km : mean24h(abs(coh_2017012)) = 0.2693 (itération1) puis 0.2489 (itération2) 0.2636
  • 2 itérations
  • METHOD1/PREPROCESS_0.1_0.2/EXP4/
2017012
FR.OGGM.00.BHZ FR.OGSM.00.BHZ FR.OGCB.00.BHN Z3.A192A.00.BHE
it1:
BF cc iter1 FR.OGGM.00.BHZ FR.OGSM.00.BHZ.png
BF cc iter1 FR.OGCB.00.BHN Z3.A192A.00.BHE.png
it2:
BF cc iter2 FR.OGGM.00.BHZ FR.OGSM.00.BHZ.png
BF cc iter2 FR.OGCB.00.BHN Z3.A192A.00.BHE.png


[edit] Distribution des δt en fonction du bin de la coherence

  • 5-10sec
  • 60min correlations
  • 2 itérations


[edit] MAX()

maxm ( corrlag5sec (GFm,lag300sec , dailystackGFm,lag300sec) )

avec m=60min

A noter que pour ces 5 jours, je n'ai trouvé aucune coherence négative dans coh(5jours,24h,14535paires,2iterations)

Binning δt = f (coherence)
2017011-2017015 2017011 2017012 2017013 2017014 2017015
iter1
BF dt distribution iter1.noabs.png
BF dt distribution 2017011 iter1 noabs.png
BF dt distribution 2017012 iter1 noabs.png
BF dt distribution 2017013 iter1 noabs.png
BF dt distribution 2017014 iter1 noabs.png
BF dt distribution 2017015 iter1 noabs.png
iter2
BF dt distribution iter2.noabs.png
BF dt distribution 2017011 iter2 noabs.png
BF dt distribution 2017012 iter2 noabs.png
BF dt distribution 2017013 iter2 noabs.png
BF dt distribution 2017014 iter2 noabs.png
BF dt distribution 2017015 iter2 noabs.png


[edit] MAX(ABS())

maxm ( ABS ( corrlag5sec (GFm,lag300sec , dailystackGFm,lag300sec) ) )

avec m=60min

Au total (120heures,14535paires,2iterations), on a 40.8% (iteration 1) et 41.5% (iteration 2) des cohérences qui sont négatives ( dont MAX(ABS()) > MAX() ).

Binning δt = f (coherence)
2017011-2017015 2017011-2017015 2017011 2017012 2017013 2017014 2017015
iter1
Coh hist iter1.abs.png
BF dt distribution iter1.abs.png
BF dt distribution 2017011 iter1.png
BF dt distribution 2017012 iter1.png
BF dt distribution 2017013 iter1.png
BF dt distribution 2017014 iter1.png
BF dt distribution 2017015 iter1.png
iter2
Coh hist iter2.abs.png
BF dt distribution iter2.abs.png
BF dt distribution 2017011 iter2.png
BF dt distribution 2017012 iter2.png
BF dt distribution 2017013 iter2.png
BF dt distribution 2017014 iter2.png
BF dt distribution 2017015 iter2.png


Si on ne sélectionne que les traces HHZ (57 traces), on a donc 57*56/2=1596 paires: (OAR_JOB_ID = 3585071:3585075)

Au total (120heures,1596paires,2iterations), on a 36.7% (iteration 1) et 36.8% (iteration 2) des cohérences qui sont négatives ( dont MAX(ABS()) > MAX() ).

Binning δt = f (coherence) -- Only HHZ pairs
2017011-2017015 2017011-2017015
iter1
Coh hist iter1.absHHZ.png
BF dt distribution iter1.absHHZ.png
iter2
Coh hist iter2.absHHZ.png
BF dt distribution iter2.absHHZ.png
[edit] MAX(ABS()) et nn_lag=150sec (au lieu de 300sec)

On pense que 300sec pour le lag de la correlation est trop élevé. Les travaux de Laurent Stehly ont montré que les niveaux de cohérence étaient bien meilleurs lorsqu'on fenétrait sur l'onde de surface.

Si l'on considère une paire relativement éloigné mais avec une bonne cohérence avec 300sec de lag:

  • Z3.A180A.00.HHZ_Z3.A206A.00.HHZ (paire indice 13398 / 14535)
  • dist=309km, azimuth 138°
  • coh moyenne sur 5days 2017011-2017015 = 0.3
  • coh median sur 5days 2017011-2017015 = 0.4
  • coh max sur 5days 2017011-2017015 = 0.75

on constate qu'il faut au moins garder un lag de 150sec.

2017012, iteration 1, Z3.A180A.00.BHZ_Z3.A206A.00.BHZ
nn_lag=300sec nn_lag=150sec
Cc iter1 Z3.A180A.00.BHZ Z3.A206A.00.BHZ 12.png
Cc iter1 Z3.A180A.00.BHZ Z3.A206A.00.BHZ 12 150.png

maxm ( ABS ( corrlag5sec (GFm,lag150sec , dailystackGFm,lag150sec) ) )

avec m=60min

Au total (120heures,14535paires,2iterations), on a 46.50% (iteration 1) et 44.73% (iteration 2) des cohérences qui sont négatives ( dont MAX(ABS()) > MAX() ).

Binning δt = f (coherence)
2017011-2017015 2017011-2017015
iter1
Coh hist iter1.abs150.png
Dt distribution iter1.abs150.png
iter2
Coh hist iter2.abs150.png
Dt distribution iter2.abs150.png

[edit] Evaluer le δt à partir du rapport des phases (alternative à la méthode xcorr)

2017012 H14 FR.OGGM.00.BHZ_FR.OGSM.00.BHZ
gauspuls et circular shift 0.9sec original data with the entire -300:+300sec GF zoom -70:50sec original data with cutted -50:+20sec GF
Gauspuls.png
FR.OGGM.00.BHZ FR.OGSM.00.BHZ.png
FR.OGGM.00.BHZ FR.OGSM.00.BHZ zoom.png
FR.OGGM.00.BHZ FR.OGSM.00.BHZ fenetrage.png
2017012 FR.OGGM.00.BHZ_FR.OGSM.00.BHZ
nn_lag method xcorr (sens-schonfelder) method phase
nn_lag=300sec
Cc iter1 FR.OGGM.00.BHZ FR.OGSM.00.BHZ 300 sens.png
Cc iter1 FR.OGGM.00.BHZ FR.OGSM.00.BHZ 300 phase.png
nn_lag=150sec
Cc iter1 FR.OGGM.00.BHZ FR.OGSM.00.BHZ 150 sens.png
Cc iter1 FR.OGGM.00.BHZ FR.OGSM.00.BHZ 150 phase.png


[edit] On fenêtre sur l'onde de surface avant la recherche des δt

Pour ce test on sélectionne reseau=FR,G,Z3 ; channel=HHZ ; lat=43-46 ; lon=3-8 ; soit 69 traces (soit 2346 paires) pour la journée 2017012.

On calcule les corrélations annuelles (=moyenne des corrélations horaires sur un an) pour ces 2346 paires.

On calcule des indices de début en fin de l'onde de surface, pour les côtés positif et négatif de la corrélation. On choisit un côté en fonction de l'enveloppe qui maximise l'amplitude.

Fig 1424.png
Fig 0038.png
Fig 1574.png
Fig 2108.png
Fig 1913.png
Fig 2588.png

La recherche de δt (avec la méthode xcorr) se fera uniquement sur cette fenêtre (de 40sec), et non plus sur les 600sec de toute la correlation.


On étudie plus en détail la paire OGGM - OGSM.

2017012 FR.OGGM.00.HHZ_FR.OGSM.00.HHZ
iter calcul δt sur les +/- 300sec de la corr calcul δt sur l'onde de surface (côté enveloppe max) calcul δt sur l'onde de surface (côté enveloppe min)
iter1
Cc iter1 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 sens.png
Cc iter1 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 SWdtP.png
Cc iter1 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 SWdtN.png
iter2
Cc iter2 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 sens.png
Cc iter2 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 SWdtP.png
Cc iter2 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 SWdtN.png
distribution des δt 2017012 FR.OGGM.00.HHZ_FR.OGSM.00.HHZ
iter calcul δt sur les +/- 300sec de la corr calcul δt sur l'onde de surface (côté enveloppe max) calcul δt sur l'onde de surface (côté enveloppe min)
iter1
Dt distribution 2017012 iter1.abs.png
Dt distribution 2017012 iter1.absSW P.png
Dt distribution 2017012 iter1.absSW N.png
iter2
Dt distribution 2017012 iter2.abs.png
Dt distribution 2017012 iter2.absSW P.png
Dt distribution 2017012 iter2.absSW N.png
DtP dtN 2017012.png
DtP dtN zoom 2017012.png
Dt over dtP 2017012.png
DtP over dtN 2017012.png

[edit] Etude statistique des cohérences négatives dans le cas MAX(ABS()) et methode xcorr

On se demande quelle est la signification de ces cohérences négatives.

Pour cette configuration m=60min, T=5-10sec, pour l'ensemble des 5 jours 2017011-2017015, à la première itération, parmi les valeurs de coh [24*5,14535], on a donc 40.8% de cohérence négative.

Proportion inv pol.png

On s'intéresse aux 6% de paires de traces dont coh < -0.4 et pour lesquelles on soupconne une inversion de polarité. On établit la liste de ces paires de traces (6% * 14535 * 120 = 103995). Et on cherche quelles paires reviennent le plus souvent dedans.

On obtient cette info :

$ cat invpol.txt
# coh  dt hour pair
-0.54  4.2  16     1
-0.43 -3.8  25     1
-0.55  4.0  27     1
-0.43  3.8  49     1
...
-0.51  2.8  65 14535
-0.47  3.0  70 14535
-0.48 -2.7  97 14535
-0.64 -3.3 103 14535
-0.43 -3.3 115 14535

On cherche combien de fois chaque paire est représentée:

$ cat invpol.txt |grep -v "#"|awk '{print $4}'|uniq -c
    11 1
    15 2
    15 3
     5 4
    17 5
...
    12 14531
    12 14532
    15 14533
    20 14534
    10 14535

Il semble que presque toutes les paires soit au moins une fois concernées... (rappel on a 14535 paires de traces différentes) :

$ cat invpol.txt |grep -v "#"|awk '{print $4}'|uniq -c|wc -l
14265

On cherche combien de paires sont concernées au moins 10% du temps (au moins 12h/120h):

$ cat invpol.txt |grep -v "#"|awk '{print $4}'|uniq -c | sort -gk 1 |awk ' $1 >= 12 '|wc -l
2404

Donc 2404/14535 = 17%

On cherche la trace qui apparaît le plus souvent:

$ for id_pair in $(cat invpol.txt |grep -v "#"|awk '{print $4}') ; do cat -n local_pair_dist.txt |grep -w "^ *${id_pair}" ; done > list_id_pair.txt
$ cat list_id_pair.txt | awk '{print $2 " " $3}'| tr -s " " "\n"|sort|uniq -c|sort -gk 1
   626 Z3.A199A.00.BHE
   670 Z3.A199A.00.BHN
   713 FR.RSL.00.BHE
   742 Z3.A199A.00.BHZ
   774 Z3.A200A.00.BHE
   790 FR.RSL.00.BHN
   798 Z3.A200A.00.BHN
   800 FR.RSL.00.BHZ
   803 FR.OGMO.00.BHN
   847 Z3.A178A.00.BHE
   853 Z3.A178A.00.BHZ
   874 Z3.A201A.00.BHN
   893 FR.BLAF.00.BHE
   897 FR.GRN.00.BHE
   906 Z3.A200A.00.BHZ
   945 Z3.A204A.00.BHZ
   948 FR.OGSM.00.BHN
   962 FR.MLYF.00.BHE
   962 FR.OGSA.00.BHE
   966 Z3.A201A.00.BHE
   970 FR.BANN.00.BHE
   984 FR.ARBF.00.BHN
   984 Z3.A181A.00.BHN
   997 Z3.A185A.00.BHE
   997 Z3.A185A.00.BHN
  1000 FR.BANN.00.BHN
  1004 Z3.A176A.00.BHE
  1008 FR.OGDI.00.BHN
  1010 FR.OGMO.00.BHE
  1012 Z3.A204A.00.BHE
  1020 FR.OGAG.00.BHN
  1022 FR.OGCB.00.BHE
  1031 FR.OGGM.00.BHE
  1032 Z3.A201A.00.BHZ
  1035 FR.OGSM.00.BHE
  1037 Z3.A188A.00.BHE
  1040 FR.OGMO.00.BHZ
  1042 FR.BLAF.00.BHN
  1057 Z3.A178A.00.BHN
  1058 FR.SAOF.00.BHE
  1062 Z3.A217A.00.BHE
  1072 FR.OGSA.00.BHZ
  1077 Z3.A176A.00.BHN
  1080 FR.OGDI.00.BHE
  1086 FR.ARBF.00.BHE
  1089 FR.SAUF.00.BHN
  1099 FR.OGSM.00.BHZ
  1101 FR.OGS3.00.BHE
  1104 FR.OGSA.00.BHN
  1114 FR.OGS2.00.BHE
  1118 FR.OGGM.00.BHN
  1123 Z3.A215A.00.BHN
  1127 Z3.A180A.00.BHE
  1129 FR.OGS2.00.BHN
  1133 Z3.A181A.00.BHZ
  1134 Z3.A188A.00.BHZ
  1138 FR.GRN.00.BHZ
  1142 FR.OGS3.00.BHN
  1142 FR.OGVG.00.BHE
  1142 Z3.A192A.00.BHE
  1146 Z3.A204A.00.BHN
  1148 FR.OGGM.00.BHZ
  1149 FR.BLAF.00.BHZ
  1152 Z3.A184A.00.BHN
  1155 FR.MLYF.00.BHZ
  1156 FR.MLYF.00.BHN
  1162 Z3.A181A.00.BHE
  1168 FR.OGAG.00.BHE
  1171 FR.TRIGF.00.BHN
  1173 FR.RUSF.06.BHE
  1173 Z3.A216A.00.BHN
  1174 FR.OGVG.00.BHN
  1176 FR.RUSF.01.BHN
  1185 Z3.A184A.00.BHE
  1186 Z3.A188A.00.BHN
  1187 Z3.A193A.00.BHE
  1191 Z3.A217A.00.BHZ
  1192 FR.RUSF.04.BHN
  1203 Z3.A193A.00.BHN
  1205 FR.FLAF.00.BHN
  1207 FR.OGS3.00.BHZ
  1207 FR.SURF.00.BHN
  1207 Z3.A206A.00.BHE
  1209 FR.OGS2.00.BHZ
  1211 FR.ESCA.01.BHE
  1212 FR.SURF.00.BHE
  1212 Z3.A215A.00.BHE
  1214 Z3.A216A.00.BHE
  1216 FR.OGDI.00.BHZ
  1216 Z3.A176A.00.BHZ
  1217 FR.TURF.00.BHN
  1222 FR.OGCB.00.BHN
  1227 FR.SAUF.00.BHE
  1239 FR.ENAUX.00.BHE
  1239 FR.ISO.00.BHE
  1240 FR.OGAG.00.BHZ
  1240 FR.OGCB.00.BHZ
  1243 Z3.A217A.00.BHN
  1248 Z3.A184A.00.BHZ
  1249 FR.ARTF.00.BHN
  1251 FR.RUSF.05.BHN
  1253 FR.SAOF.00.BHN
  1254 FR.SAOF.00.BHZ
  1255 FR.GRN.00.BHN
  1255 FR.RUSF.05.BHZ
  1256 Z3.A193A.00.BHZ
  1258 Z3.A180A.00.BHN
  1259 FR.BANN.00.BHZ
  1263 FR.RUSF.03.BHE
  1266 FR.ESCA.01.BHN
  1269 FR.RUSF.03.BHZ
  1271 FR.SPIF.00.BHN
  1273 Z3.A205A.00.BHN
  1273 Z3.A206A.00.BHN
  1274 Z3.A185A.00.BHZ
  1276 FR.RUSF.04.BHZ
  1279 FR.RUSF.06.BHZ
  1283 FR.RUSF.07.BHN
  1287 FR.RUSF.07.BHE
  1287 FR.RUSF.07.BHZ
  1291 FR.SPIF.00.BHZ
  1292 FR.RUSF.01.BHZ
  1296 FR.RUSF.06.BHN
  1304 FR.OGVG.00.BHZ
  1312 FR.CALF.00.BHN
  1323 FR.TRIGF.00.BHE
  1325 FR.ISO.00.BHN
  1329 FR.RUSF.04.BHE
  1336 FR.OGS1.00.BHZ
  1339 FR.RUSF.03.BHN
  1339 FR.SURF.00.BHZ
  1345 Z3.A180A.00.BHZ
  1346 Z3.A206A.00.BHZ
  1346 Z3.A215A.00.BHZ
  1347 FR.ARBF.00.BHZ
  1353 FR.RUSF.01.BHE
  1355 Z3.A194A.00.BHN
  1361 FR.FLAF.00.BHZ
  1365 Z3.A192A.00.BHN
  1368 Z3.A216A.00.BHZ
  1371 FR.SAUF.00.BHZ
  1375 Z3.A194A.00.BHE
  1383 FR.RUSF.05.BHE
  1386 FR.ARTF.00.BHE
  1387 FR.CALF.00.BHE
  1389 FR.SPIF.00.BHE
  1392 FR.ISO.00.BHZ
  1400 FR.FLAF.00.BHE
  1418 FR.CALF.00.BHZ
  1420 FR.ESCA.01.BHZ
  1443 FR.ENAUX.00.BHN
  1443 FR.MVIF.00.BHN
  1448 FR.ENAUX.00.BHZ
  1465 FR.EILF.00.BHN
  1465 Z3.A205A.00.BHZ
  1471 FR.BSTF.00.BHE
  1485 FR.MVIF.00.BHE
  1486 FR.ARTF.00.BHZ
  1488 FR.OGS1.00.BHE
  1511 FR.TRIGF.00.BHZ
  1518 FR.EILF.00.BHE
  1566 FR.OGS1.00.BHN
  1635 FR.MVIF.00.BHZ
  1640 FR.EILF.00.BHZ
  1653 FR.BSTF.00.BHN
  1668 FR.BSTF.00.BHZ
  1682 Z3.A194A.00.BHZ
  1733 Z3.A192A.00.BHZ
  1929 Z3.A205A.00.BHE
  2245 FR.TURF.00.BHE
  2833 FR.TURF.00.BHZ

On obtient donc en fin de liste les traces avec un indice qui indique le nombre de fois ou cette trace est apparue dans paire de trace considérée comme suspecte sur une heure donnée entre 2017011 et 2017015.


On cherche LA ou LES paire la plus représentée:

$ cat invpol.txt |grep -v "#"|awk '{print $4}'|uniq -c | sort -gk 1
     1 10000
     1 10006
     1 10012
     1 10031
     1 10046
...
    37 12198
    37 2205
    37 8037
    38 12523
    39 12825
    41 12331
    42 7808
    44 4038
    44 4182
    44 9325
$ cat -n local_pair_dist.txt |grep -w 4038
 4038	FR.ESCA.01.BHN FR.TURF.00.BHZ 16182.162
$ cat -n local_pair_dist.txt |grep -w 4182
 4182	FR.ESCA.01.BHZ FR.TURF.00.BHZ 16182.162
$ cat -n local_pair_dist.txt |grep -w 9325
 9325	FR.OGSA.00.BHZ FR.TURF.00.BHE 142000.797
$ 

Donc on a affaire à 3 paires de traces distantes de 16km et 142km. On regarde le détail du calcul aux 2 premières itérations pour ces paires de traces (EXP5):

$ h5ls /scratch_md1200/lecointre/IWORMS/METHOD1/PREPROCESS_0.1_0.2/EXP5/2017.011.trace.int8.h5
FR.ESCA.01.BHN           Dataset {864000}
FR.ESCA.01.BHZ           Dataset {864000}
FR.OGSA.00.BHZ           Dataset {864000}
FR.TURF.00.BHE           Dataset {864000}
FR.TURF.00.BHZ           Dataset {864000}
  • m=60' : OAR_JOB_ID: 3542937 ... 3542941 (3 itérations pour voir les GFs et la GFstack après 2 corrections successives)
  • m=120' : OAR_JOB_ID: 3542944 ... 3542948
FR.ESCA.01.BHN_FR.TURF.00.BHZ
iteration 1 iteration 2 iteration 3 ( sum(abs(coh(iter1:iter2))) et sum(dt(iter1:iter2)) )
2017011 (m=60)
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ 11.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter2 11.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter3 11.png
2017012 (m=60)
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ 12.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter2 12.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter3 12.png
2017013 (m=60)
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ 13.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter2 13.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter3 13.png
2017014 (m=60)
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ 14.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter2 14.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter3 14.png
2017015 (m=60)
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ 15.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter2 15.png
Invpol EXP5 FR.ESCA.01.BHN FR.TURF.00.BHZ iter3 15.png
FR.ESCA.01.BHZ_FR.TURF.00.BHZ
iteration 1 iteration 2 iteration 3 ( sum(abs(coh(iter1:iter2))) et sum(dt(iter1:iter2)) )
2017011 (m=60)
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ 11.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 11.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 11.png
2017011 (m=120)
Invpol EXP5 120 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter1 11.png
Invpol EXP5 120 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 11.png
Invpol EXP5 120 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 11.png
2017012 (m=60)
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ 12.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 12.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 12.png
2017012 (m=120)
Invpol EXP5 120 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter1 12.png
Invpol EXP5 120 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 12.png
Invpol EXP5 120 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 12.png
2017013 (m=60)
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ 13.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 13.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 13.png
2017014 (m=60)
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ 14.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 14.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 14.png
2017015 (m=60)
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ 15.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter2 15.png
Invpol EXP5 FR.ESCA.01.BHZ FR.TURF.00.BHZ iter3 15.png
FR.OGSA.00.BHZ_FR.TURF.00.BHE
iteration 1 iteration 2 iteration 3 ( sum(abs(coh(iter1:iter2))) et sum(dt(iter1:iter2)) )
2017011 (m=60)
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE 11.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter2 11.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter3 11.png
2017012 (m=60)
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE 12.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter2 12.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter3 12.png
2017013 (m=60)
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE 13.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter2 13.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter3 13.png
2017014 (m=60)
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE 14.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter2 14.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter3 14.png
2017015 (m=60)
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE 15.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter2 15.png
Invpol EXP5 FR.OGSA.00.BHZ FR.TURF.00.BHE iter3 15.png
  • Robustesse 60' / 120' ? ...
  • On fait parfois apparaître un δt qui n'était pas visible au départ ...
  • Le +/-δt=+/-T/2 induit par une cohérence très négative n'est pas toujours bien compensé par un -/+δt=-/+T/2 à l'itération suivante ...

[edit] Masquer les Δt[traces_id=171] finaux à partir des cohérences des paires (coh[2,paires_id=14535]) de traces associées ?

  • Calcul complet jusqu'au Δt (avec méthode C.Picard, donc pas de seuil sur cohérence: on garde tout pour résoudre le système linéaire)
  • On recherche les cohérence max avec MAX(ABS())
  • T=5-10sec
  • Deux expériences:
    • 5jours, m=60' -> 120 timesteps
    • 5jours, m=120' -> 60 timesteps

Pb1 : Pour chaque trace_id, on a peu de paires associées (=parmi les 170 voisins de cette trace_id) qui respectent le critère de cohérence (meanit1,it2(abs(coh))>=seuil). On a testé seuil=0.3, 0.4, 0.5.

Pb2 : Ce masque semble assez variable en temps...

2017011-2017015
mean(abs(coh))>=0.3 mean(abs(coh))>=0.4 mean(abs(coh))>=0.5
m=60' -> 120 timesteps
Rate voisins seuilcoh 0.3.png
Rate voisins seuilcoh 0.4.png
Rate voisins seuilcoh 0.5.png
m=120' -> 60 timesteps
Rate voisins seuilcoh 0.3 2h.png
Rate voisins seuilcoh 0.4 2h.png
Rate voisins seuilcoh 0.5 2h.png

On fixe:

  • 5jours, m=120' -> 60 timesteps
  • cohseuil=0.4
  • >=50% de paires associées qui valident le critère de seuil de cohérence
méthode AtA=NId puis masque là ou plus de 50% de voisins avait une coh < seuil on supprime des lignes de A et on inverse A (matlab pinv)
Dt 6.png
Dt inv 6.png
Dt2 inv 6.png
Dt 5.png
Dt inv 5.png
Dt2 inv 5.png
Dt 4.png
Dt inv 4.png
Dt2 inv 4.png
Dt 3.png
Dt inv 3.png
Dt2 inv 3.png
Dt 2.png
Dt inv 2.png
Dt2 inv 2.png
Dt 1.png
Dt inv 1.png
Dt2 inv 1.png
Slmask.png
Slmask.png

[edit] Correlations

[edit] Method

  • fréquence: 30min
  • toutes les paires, quelque soit l'interdistance: 171 traces -> N(N-1)/2 pas d'autocorr -> 14,535 correlations/30min
  • lag = 300sec (ToDo: Lag = f (interdistance stations) = ( distance / vitesse ) x 1.5 avec vitesse=2km/s)
  • save as HDF5 f32 compressed dataset (gzip or scale offset compression ?)
  • corr_bypair séquentiel ou corr_planmany avec FFTW multithreadée si mémory needs ? Le projet SanJacinto a montré de meilleures perfs avec l'option _planmany (quitte à soumettre sur plusierus coeurs quand les besoins mémoire le nécessitent et à activer la FFTW multithreadée), mais la version _bypair est plus simple (code, déploiement, ...). Pour la suite des calculs (stack des Green Functions pour toutes les paires), le code planmany sera plus approprié (avec boucle la plus externe sur le temps : par pas de 30min).

[edit] Source code

$ svn co https://forge.osug.fr/svn/isterre-correlation/IWORMS/IWORMS_correlation_src IWORMS_correlation_src

A    IWORMS_correlation_src/submit_grid_bypair.sh
A    IWORMS_correlation_src/corr_bypair.f90
A    IWORMS_correlation_src/corr_planmany.f90
A    IWORMS_correlation_src/namelist.skel
A    IWORMS_correlation_src/Makefile
Révision 212 extraite.
$

[edit] Compilation

$ module load intel-devel/13 hdf5/1.8.14_intel-13.0.1 zlib/1.2.7_gcc-4.6.2 szip/2.1_gcc-4.6.2 fftw/3.3.4-precise_intel-13.0.1
$ make corr_bypair
h5pfc -cpp -Dlk_scaleoffset -Dlk_float -O2  -o corr_bypair corr_bypair.f90 -I/applis/ciment/v2/stow/x86_64/intel_13.0.1/fftw_3.3.4-precise/include /applis/ciment/v2/stow/x86_64/intel_13.0.1/fftw_3.3.4-precise/lib/libfftw3f.a /applis/ciment/v2/stow/x86_64/intel_13.0.1/fftw_3.3.4-precise/lib/libfftw3f_threads.a
$ make corr_planmany
h5pfc -cpp -Dlk_scaleoffset -Dlk_float -O2  -o corr_bypair corr_bypair.f90 -I/applis/ciment/v2/stow/x86_64/intel_13.0.1/fftw_3.3.4-precise/include /applis/ciment/v2/stow/x86_64/intel_13.0.1/fftw_3.3.4-precise/lib/libfftw3f.a /applis/ciment/v2/stow/x86_64/intel_13.0.1/fftw_3.3.4-precise/lib/libfftw3f_threads.a

[edit] Computing

Namelist

Submission:

$ oarsub -S "./submit_grid_bypair.sh 2017 011"

Performance:

Performance calculs correlations corr_bypair vs corr_planmany
corr_bypair corr_planmany
node luke4 (1core, seq) luke4 (1core, seq)
IO (read and write) /nfs_scratch /nfs_scratch
prec single precision computation: compile with activate lk_float CPPkey single precision computation: compile with activate lk_float CPPkey
output F32 with lk_scaleoffset CPPkey F32 (h5 lite interface)
memory 38MB 470MB
CPUcost 17 CPU-min (0.070sec/1day-pair -> 0.0015sec/30min/pair) : this implies a writing disk operation every 0.070sec -> problem if deploying on cigri. but we don't want to deploy corr code only and write all correlations, isn't it ? 6 CPU-min (7.26sec/30min-allpairs -> 0.0005sec/30min/pair) : this implies a writing disk operation every 0.0005sec -> problem if deploying on cigri. but we don't want to deploy corr code only and write all correlations, isn't it ?

Outputs:

$ ls -lh 2017011.000000.h5 
-rw-r--r-- 1 lecointre l-isterre 7.5G Jun  7 20:58 2017011.000000.h5
$ h5ls -v 2017011.000000.h5
Opened "2017011.000000.h5" with sec2 driver.
FR.ARBF.00.BHE_FR.ARBF.00.BHN Dataset {48/48, 6001/6001}
    Location:  1:800
    Links:     1
    Chunks:    {48, 6001} 1152192 bytes
    Storage:   1152192 logical bytes, 576118 allocated bytes, 199.99% utilization
    Filter-0:  scaleoffset-6 OPT {0, 5, 288048, 1, 4, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    Type:      native float
FR.ARBF.00.BHE_FR.ARBF.00.BHZ Dataset {48/48, 6001/6001}
    Location:  1:580142
    Links:     1
    Chunks:    {48, 6001} 1152192 bytes
    Storage:   1152192 logical bytes, 540112 allocated bytes, 213.32% utilization
    Filter-0:  scaleoffset-6 OPT {0, 5, 288048, 1, 4, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    Type:      native float
FR.ARBF.00.BHE_FR.ARTF.00.BHE Dataset {48/48, 6001/6001}
    Location:  1:1123150
    Links:     1
    Chunks:    {48, 6001} 1152192 bytes
    Storage:   1152192 logical bytes, 540112 allocated bytes, 213.32% utilization
    Filter-0:  scaleoffset-6 OPT {0, 5, 288048, 1, 4, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
    Type:      native float
...

[edit] Calcul des δt

Loop on days --> Cigri parallelisation
  Loop on pair of traces
    Si dist < nn_dist
      Computes houry-correlations with 300s lag -> 24 1h-GF : {24,6001}
      Loop on iterations (nn_iter=2)
        Computes 1d-stackedGF, starting from 1h-GFs -> {6001}
        Loop on hours (nn_hour=24)
          Compute corr(1d_GF,1h_GF): correlations between 1day-stacked_GF and the 1h-GF with 5s lag --> {101}
          Si MAX(ABS(corr(1d_GF,1h_GF))) > nn_coherence_seuil
            δtiter_n = MAXLOC(ABS(corr(1d_GF,1h_GF)))    ! ABS ou pas ???
            Update 1h-GF (circular shift) 
          Sinon
            ?
          Fin Si nn_coherence_seuil
        End Loop on hours
      End Loop on iterations
      δt = δtiter_1 + δtiter_2 --> {24}
    Fin Si dist 
  End Loop on pair of traces
End Loop on days

[edit] Calcul des Δt et Tests avec les profils PE/PF du réseau X7 (Pyrope)

Test méthode AT*A=N*Id avec hypothèse de fermeture (Σ Δt = 0) avec dataset Pyrope.

On dispose de 2941 δt obtenus à partir d'un dataset de 103 stations : 69 Pyrope (PY) et 34 profils PE/PF.

En fait on ne dispose que des paires PE/PF avec PY ainsi que des paires PE/PF avec PE/PF (y compris les autocorr), mais on n'a pas PY avec PY :

  • ( corr PE/PF - PY ) + ( corr PE/PF - PE/PF ) = ( 34 * 69 ) + ( 34*35/2 ) = 2941

On sélectionne uniquement les 34 traces PE/PF pour lesquelles on a toutes les paires, sans les autocorr. On a donc 33*34/2 = 561 paires.

On construit AT[N,N(N-1)/2+1] avec N=34 et en respectant l'hypothèse de fermeture.

On calcule, pour chaque trace js, Δt(js) = (1/N) * < AT(js,:),δt(:) >

Δt sur profil PE/PF
1day 5days
03-05sec
Dt pyrope.db dt stack 01 03-05s win type sw.mat.png
Dt pyrope.db dt stack 05 03-05s win type sw.mat.png
03-10sec
Dt pyrope.db dt stack 01 03-10s win type sw.mat.png
Dt pyrope.db dt stack 05 03-10s win type sw.mat.png
05-07sec
Dt pyrope.db dt stack 01 05-07s win type sw.mat.png
Dt pyrope.db dt stack 05 05-07s win type sw.mat.png
07-10sec
Dt pyrope.db dt stack 01 07-10s win type sw.mat.png
Dt pyrope.db dt stack 05 07-10s win type sw.mat.png
10-20sec
Dt pyrope.db dt stack 01 10-20s win type sw.mat.png
Dt pyrope.db dt stack 05 10-20s win type sw.mat.png

On se focalise sur la bande de fréquence 3-5sec, pour les deux expériences 1day et 5days.

Δt 03-05sec sur profil PE/PF
1day 5days
Σ δt ≠ 0 (comme ci dessus)
Dt pyrope.db dt stack 01 03-05s win type sw.mat.png
Dt pyrope.db dt stack 05 03-05s win type sw.mat.png
Σ δt
Dt sum.db dt stack 01 03-05s win type sw.mat.png
Dt sum.db dt stack 05 03-05s win type sw.mat.png
On force Σ δt = 0 (on applique ∀ (i,j) ∈ {Npair} , δti,j=δti,j-Σ δt / Npair)
Dt2 pyrope.db dt stack 01 03-05s win type sw.mat.png
Dt2 pyrope.db dt stack 05 03-05s win type sw.mat.png
2D plots ( Σ δt ≠ 0 / Σ δt = 0)
Dt 2Dplot 1.db dt stack 01 03-05s win type sw.mat.png
Dt 2Dplot 1.db dt stack 05 03-05s win type sw.mat.png
2Dplots (suite)
Dt 2Dplot 2.db dt stack 01 03-05s win type sw.mat.png
Dt 2Dplot 2.db dt stack 05 03-05s win type sw.mat.png
2Dplots (suite)
Dt 2Dplot 3.db dt stack 01 03-05s win type sw.mat.png
Dt 2Dplot 3.db dt stack 05 03-05s win type sw.mat.png
2Dplots (suite)
Dt 2Dplot 4.db dt stack 01 03-05s win type sw.mat.png
Dt 2Dplot 4.db dt stack 05 03-05s win type sw.mat.png
2Dplots (suite)
Dt 2Dplot 5.db dt stack 01 03-05s win type sw.mat.png
Dt 2Dplot 5.db dt stack 05 03-05s win type sw.mat.png
coeff associés aux δt
Coeffdt distribution.db dt stack 01 03-05s win type sw.mat.png
Coeffdt distribution.db dt stack 05 03-05s win type sw.mat.png
coeff associés aux δt
Histcoeffdt distribution.db dt stack 01 03-05s win type sw.mat.png
Histcoeffdt distribution.db dt stack 05 03-05s win type sw.mat.png

Il y a ~50,000 valeurs de coeff (sur ~300*561=170,000) qui sont nulles. Elles correspondent toutes à des δt=0.

Les sauts brutaux de Δt correspondent aux dates suivantes:

  • 2011-12-15 : PF06, PF09, PF15
  • 2012-01-24 : PF09, PF15
  • 2012-02-24 : PF01, PF12, PF15
  • 2012-04-19 : PF01, PF03, PF09
  • 2012-05-11 : PF01, PF09, PF12
  • 2012-07-02 : toutes
  • 2012-08-02 : PF03

Il semble que ces sauts soient tout aussi brutaux dans le cas du stack 5days que dans le cas du stack 1day.

[edit] Compte rendu de réunions

[edit] 05/04/2017 : Réunion de lancement du projet

Participants: P. Roux, C. Pequegnat, A. Lecointre, V. Louvet, B. Bzeznik, P. Boué, ... ?

AO CNRS ne finance que 6000euros sur les 22000 demandés. possibilité de mutualiser avec les 20kE d'Imag'in renouvelés en 2017 (qui avaient permis d'acquérir 2 noeuds luke24 et luke25 en 2016)

Nouveaux participants potentiels au projet : O. Henrot, P.-A. Bouttier

Quel dataset ?

  • AlpArray ?
  • VolcArray ?
  • Garambois/Voisin ?

Corrélations ?

  • journalières ?
  • horaires ?

1 paire de stations -> 9 corrélations (car 3 composantes) -> Moyenne des 9 ou on garde les 9 pour le calcul des δt ?

Preprocess : Plateforme BigData GRICAD ?

Flux continu ? utiliser les données de /media/resif (données corrigées) à resif9@BIO au SIMSU (actualisées 2 fois par jour)

Comment prendre en compte l'ajout complémentaire progressif de données ?

Schéma complet du workflow : ToDo (-> Albanne)

Besoin de demande de conseils aux collègues de MathAppli pour l'inversion de la matrice G. Violaine prend contact et programme une réunion.

[edit] 09/05/2017 : Réunion Choix des données et prétraitements

Participants : P. Roux, F. Brenguier, C. Pequegnat, L. Stehly, A. Lecointre

Réseau:

  • AlpArray: 121 stations -> on en choisit 50 sur critère d'interdistance minimum : zone de densité max: Alpes maritimes
  • 3 composantes -> jeu test d'~150 traces
  • ~100Hz : grep station code = HH[ENZ]

Catherine explique à Albanne le fonctionnement du portail web RESIF pour la récupération du fichier ascii stationsByCodes. Pb de retour à la ligne:

sed 's/^M/\n/g'

-> Catherine fait un ticket au RESIFportal admin.

Preprocess: discussion autour de deux méthodes à tester...

[edit] 15/05/2017 : Réunion Prétraitement des données

Participants: G. Cougoulat, O. Henriot, C. Lenne, V. Louvet, B. Bzeznik, A. Lecointre, ... ?

Les données AlpArray nécessitent d'être nettoyées avant le preprocess:

  • Il n'y a pas que du 100Hz : on trouve aussi du 125Hz, du 200Hz, ...
  • Il faut recaler sur des centièmes de secondes pleines. Ce décalage semble fixe pour un segment de trace donné. mais pas forcément pour deux segments consécutifs (voir AJAC entre 2015 et 2016...). Et ce décalage n'est pas forcément fixe pour les trois canaux d'une même station (voir ARTF)
  • Ca déborde sur J+1. Est il possible que ca déborde aussi sur J-1 ?
  • Il faut parfois aller chercher le début de J dans J-1 ? on ne peut pas travailler indépendamment d'une journée à l'autre. il faudra télécharger 3 jours pour pretraiter un jour (IOread *3)
  • Peut on avoir des recouvrement non égaux d'un jour à l'autre : question non encore élucidée ?

Il faudra bien prendre en compte le network code ainsi que le location code (en plus du station code et channel code) pour l'identification des traces (4 identifiants pour une trace unique).

Glenn pose la question de l'éventuelle nécessité de déconvoluer (enlever la réponse instrumentale) avant tout travail avec les données. Finalement il semble que ce ne soit pas nécessaire: [...] techniquement si tu regardes des variations de temps au cours du temps par couple tu ne devrais pas avoir besoin de deconvoluer. [...], je regarde pour te donner tous les scripts que j'ai et surtout voir si par requete sur webservice on s'affranchira pas des problèmes de recalage des journées... en requetant 1an directement par station et découper par obspy en plus c'est plus malin pour le map reduce. A voir [...] -> en attente de complément de réponse.

Violaine souligne l'intérêt de développer une approche comparative Grille / BigData. Albanne développerait les outils à déployer sur la grille et GRICAD fournirait un ingénieur du LJK pour développer tout ou partie des outils à déployer sur une plateforme BigData avec MapReduce ( quel FS distribué ? : HDFS Hadoop ?, Spark ?, ... ?).

[edit] 18/05/2017 : Réunion Inversion matrice

Participants: P. Roux, L. Stehly, A. Lecointre (à distance), V. Louvet, C. Picard, ... ?

Présentation slides Laurent Stehly : demander autorisation pour intégrer slides à ce wiki ?

Objectif scientifique: Soit un réseau avec n=1000 stations, évaluer si les 1000 stations sont bien synchronisées en temps ou non ? Si non, être capable de calculer un Δt pour chaque n et corriger (=shifter) le défaut horaire sur les données.

Inversion matrice [n,n(n-1)/2]:

  • 6MB (r4) si n=150 et si on stocke tout
  • mais matrice creuse: 2 entrées non nulles par ligne
  • méthode directe ?
  • approche par blocs ?

Calcul des corrélations :

  • tri de stations selon ordre spatial
  • chercher un jour avec au moins 150 traces valides
  • quelle est la bonne échelle de temps ? 30min ?
  • quel lag ? lapse time pour le calcul de la fonction de Green = 2 * distance / V avec V=2km/s

Choisir selon critère les stations qui vont intervenir dans l'inversion (critère de fiabilité), pouvoir enlever des corrélations, ne stacker que les corr qui se ressemblent.

Laurent Stehly: Si l'on fait l'hypothèse qu'un certain nombre de stations ne dérivent pas (raisonnable dans le cas d'un réseau permanent, déraisonnable dans le cas de stations OBS Ocean Bottom Seismometer), on peut déduire l'erreur d'horloge directement en comparant la médiane globale d'une station par rapport à toutes les autres, pas d'inversion.

[edit] 20/06/2017 : Réunion Inversion matrice

Participants: P. Roux, L. Stehly, A. Lecointre, C. Picard, L. Métivier

1) Il serait possible de construire une solution sans résoudre le système linéaire.

On a un graphe orienté.

La matrice de passage (notée G dans l'article Schonfelder). On la note A. Exemple si n=4 traces (4 stations avec une seule composante):

     1 -1  0  0
     1  0 -1  0
     1  0  0 -1
A =  0  1 -1  0
     0  1  0 -1
     0  0  1 -1 

ATA = une matrice [n,n] avec n-1 sur la diagonale et -1 partout ailleurs. Cette matrice de laplacien, on connaît ses valeurs et vecteurs propres:

  • λ = 0 : de multiplicité 1
  • λ = n : de mult n-1 -> EV de dim n-1

On peut donc diagonaliser ATA sous la forme PDP-1

2) Hypothèse supplémentaire: relation de fermeture Σ δ ti,j = 0

Cela revient à ajouter une dernière ligne de 1 à la matrice A. Dans ce cas ATA devient diagonale = n*Id

Quand on a un graphe orienté complet, la matrice est donc trivialement inversible.

3) Remarques:

  • Il est impossible de quantifier l'erreur avec cette méthode -> il faudrait avoir ATCA avec C matrice de covariance d'erreur d'observation (barre d'erreur sur chaque trajet) qui serait liée à la cohérence (= coeff de ressemblance) entre les 2 fonctions de Green dailystack et horaire.
  • Il est impossible de n'avoir que 2 traces

4) En réalité on n'a jamais un ensemble complet. on aura toujours un certain nombre de δt entachés d'erreur. Ce ne sont pas des zéros, ce sont des valeurs indéterminées (car le coefficient de ressemblance < seuil).

  • analyser le graphe (pour chaque heure !) pour en trouver des sous ensembles complets ? et de proche en proche recontruire la solution complète ? Comment identifier un sous graphe complet ? (-> voir si question abordée en Théorie des graphes)
  • dégressif: d'abord dailyGF versus yearlyGF, puis hourlyGF versus dailyGF ?

5) Exemple avec n=4 et δt1,3 indéterminé et relation de fermeture, alors ATA =

       3 0 1 0
A'A =  0 3 0 0
       1 0 3 0
       0 0 0 4
  • N'est elle pas à nouveau trivialement inversible ? Puisque chacun des Det d'ordre n-1 se retrouve diagonal à son tour ?
  • Que se passe t il quand plusieurs δti,j sont indéterminés ?
  • Que se passe t il quand on a un i tel que les δti,k sont tous indéterminés pour tout k ? -> dans ce cas A devient singulière donc on supprime totalement une colonne de A et on abaisse l'ordre du système linéaire à résoudre ? On n'aura donc jamais le résultat de Δti. Par contre on retrouve ATA = (n-2)*Id

6) Laurent fournit à Albanne un jeu de données δti,j Pyrope pour lequel on connaît à l'avance Δti. Pour tester l'approche 2)

[edit] 21/07/2017 : Réunion Qualité des Données et Webservices

Participants : V. Louvet, G. Cougoulat, A. Lecointre

1) Point méthodologie / qualité des données:

Albanne expose les difficultés de définir une méthodologie de référence pour le calcul des δt :

  • pic de δt autour de T/2 : lié au problème des "inversions de polarité" ?
  • valeurs louches de δt pour les 1 ou 2 premières heures de la journée : à confirmer ?

Glenn propose d'éliminer les stations de mauvaise qualité pour définir la méthodologie sur des données "propres". Il s'agit en gros de ne retenir que le réseau Z3 ainsi que les stations "de qualité Z3" dans les réseau FR/RD. Ceci impose de se déplacer en août 2016 pour effectuer les tests car pas de Z3 en 2017 ? Pourquoi ai je noté cela ?. Glenn fournit une liste des "bonnes stations" FR/RD/G : File:Liste stations permanentes.pdf. On conserve les appareils: Taurus, Q330, Q330HR, Q330S, CEA. On supprime les appareils: Kephren, Agecodagis_Titan, Staneo_D6BB-DIN, CMG_DM24, Titan3.

$ cat -n liste_stations_FR_RD_G_OK.txt
    1	FR BANN
    2	FR BOUC
    3	FR CAMF
    4	FR CFF
    5	FR CHIF
    6	FR CHMF
    7	FR CLAF
    8	FR CLEV
    9	G CLF
   10	FR DOU
   11	G ECH
   12	FR ENAUX
   13	FR FLAF
   14	FR GRN
   15	RD LOR
   16	RD MFF
   17	FR MOF
   18	FR MORSI
   19	RD MTLF
   20	FR OG02
   21	FR OGCB
   22	FR OGMO
   23	FR OGMY
   24	FR OGRV
   25	FR OGS1
   26	FR OGSA
   27	FR OGSI
   28	FR OGSM
   29	FR OGVG
   30	RD ORIF
   31	RD PGF
   32	FR RIVEL
   33	FR RONF
   34	RD ROSF
   35	G SSB
   36	FR STR
   37	FR TRBF
   38	FR TRIGF
   39	FR WLS
$ cat -n liste_stations_FR_RD_G_NOK.txt
    1	FR AJAC
    2	FR ARBF
    3	FR ARTF
    4	FR ASEAF
    5	FR ATE
    6	FR BLAF
    7	FR BSTF
    8	FR CALF
    9	FR CORF
   10	FR EILF
   11	FR ESCA
   12	FR FILF
   13	FR FNEB
   14	FR ISO
   15	FR LRVF
   16	FR MLS
   17	FR MLYF
   18	FR MON
   19	FR MONQ
   20	FR MVIF
   21	FR OG35
   22	FR OGAG
   23	FR OGDI
   24	FR OGGM
   25	FR PAND
   26	FR PYLO
   27	FR RENF
   28	FR RSL
   29	FR RUSF
   30	FR SAOF
   31	FR SAUF
   32	FR SJAF
   33	FR SMPL
   34	FR SPIF
   35	FR SURF
   36	FR TERF
   37	FR TURF
   38	FR UNCO
   39	FR URDF
   40	FR VIEF
$ 

2) Point approche BigData pour preprocess:

Violaine rappelle l'intérêt de tester l'approche comparative Grille / BigData pour le prétraitement des données. Glenn se fait ouvrir un compte au LAL pour l'utilisation d'une plateforme Spark/HDFS.

Le code de prétraitement n'est actuellement capable que de traiter des données présentes sur le NFS /resif (accessible depuis luke4 et luke24). Albanne adapte le code pour offrir la possibilité de requeter des données par webservices ws.resif.fr: les 2 options suivantes doivent être ajoutées:

  • wget GET method
  • obspy client fdsn

-> Fait le 27/07/2017

Cette requete par webservice ouvre le code à une utilisation plus large que RESIF.

On se demande si la requete par webservice plutôt que la lecture de fichiers depuis NFS /resif résoudra les problèmes de recalage du vecteur temps. -> la mise en oeuvre au 27/07/2017 montre que ce n'est pas le cas.

On se demande si la requete par WS peut permettre de s'affranchir totalement de la notion de journée: Appliquer le prétraitement directement sur 1h, 1mois, 1an ? mais comment gérer le gapfill dans ces conditions ?

[edit] 28/08/2017 : Réunion Faisabilité approche BigData pour preprocess

Participants : G. Cougoulat, C. Cancé (gr GRICAD Pôle de données), A. Lecointre

Durée : 1h15

  • Thématique sismologie, équipe ondes.
  • Problématique : Souhait équipe d'une approche comparative pour du prétraitement de données: hpc vs hadoop map-reduce
  • situation actuelle:
    • Analyse périodique de 150 traces journalieres issues de capteurs sur 3 composantes chacun (au maximum).
    • usage de luke4 ou luke24
    • traitement en 2 phases
1. pré-traitement de ces vecteurs (python)
-recalage de l'origine
-recalage frequence échantillonage (100hz) puis decimation
taille vecteur inital si 100 hz: 3600*100*24= 8 640 000 entrées / jour (valeurs en double précision float64) => 69 120 000 octets => 70 Mo (fichiers de ~30MB en steim). De plus après décimation on travaille sur du 10Hz (vecteurs de 7Mo).
temps d'execution (io incluses): de l'ordre de 5 sec (usage de WS echanges avec Resif ou depuis filesystem (6s))
=> traitement aboutit à un vecteur de valeurs au choix {+1,-1}.
actuellement consommation RAM du traitement peut atteindre 2 Go (à creuser ?)
2. phase de correlations entre les vecteurs des composantes enregistrées (n*(n-1)/2) calculs : 15min actuellement
  • Situation envisagée à terme:
    • augmentation du nombre de capteurs à plusieurs milliers
    • augmentation de la fréquence d'échantillonage ( > à 250 hz ds le projet resolve )

=> quels moyens seront à terme les mieux adaptés ces traitement de calculs à priori simples sur de "grands" vecteurs de données ?

  • Propositions actions:
    • 1) exploration des moyens de traitement en priorité "in memory"
- plateforme "LAL" in2p3 infrastructure nationale propose environnement spark hdfs avec python
- autre piste de "cluster in memory" type elasticsearch à creuser (ex:https://www.elastic.co/fr/blog/elasticsearch-as-a-time-series-data-store), peut être plus adapté à notre problématique actuelle "middle big data"
=> évaluation de la limite
    • 2) exploration de hadoop ds un second temps ?
  • pas d'urgence particulière -> notre cas sera évoqué en réunion Pôle de Données jeudi 31 aoüt matin.


[edit] 06/09/2017 : Réunion Inversion de Polarité , Etude des Figures , Ebauche de stratégie pour seuil cohérence

Participants : P. Roux, A. Lecointre

Durée : 1h30

  • Il y a eu une erreur dans certaines expériences / figures du wiki. Les Green Functions ont toujours été calculées sur des segments de durée 30min, quelquesoit le paramètre m = 30min , 60min , 120min. On a utilisé uniquement les 30 premières minutes de données (nn_unitseg=1800sec n'a pas été updaté dans la namelist), bien qu'on ait bien bouclé par pas de 30min , 60min ou 120min. Bug corrigé, expériences refaites, les figures ont été mises à jour. Cela modifie peu les figures. On note un léger élargissement vers la droite (vers les cohérences plus élevées) des cloches de distributions des cohérences.
  • Figures "Cohérence des corrélations sur la journée":
    • il faut enlever la moyenne sur 5 jours et faire apparaître tous les points
    • il faut faire apparître distinctement les coh négatives et positives
    • Objet MatLab h=histogram() pour récupérer les valeurs pour chaque bin et faire une représentation sous forme de courbe (permettra de superposer plusieurs représentations: par exemple les 3 bandes de fréquences).
    • On dirait que nos distributions suivent une loi de Poisson. A vérifier (MatLab histfit avec paramètre dist='poisson'). Confère GJI2013 Fig4 : Monitoring fault zone environments with correlations of earthquake waveforms File:GJI-2013-Roux-gji ggt441.pdf. Validerait l'aspect aléatoire entre 2 évènements successifs (ce qu'on attend).
  • Figures "Distribution des δt en fonction du bin de la coherence" :
    • Tant qu'à refaire des figures, on améliore certaines représentations : on distingue les cohérences positives ou négatives dans le cas du MAX(ABS()). Avant on plottait des valeurs absolues et tout était mélangé. On remarque des grappes de δt à +/- T/2 qui correspondent à des cohérences négatives : cad avec MAX(ABS()) > MAX(). Qu'est ce que c'est ?
    • Ca concerne 40% des paires de traces !
    • En réalité certaines (beaucoup ?) de ces δt à +/- T/2 sont détectées pour un niveau de cohérence faible.
    • On pickera quelques paires de traces concernées (avec forte coherence = -0.8). Plotter le détail du suivi de ces expériences sur 120h.
  • Pb des autocorr:
    • ne pas prendre en compte les autocorrélations (comme on le fait actuellement) : pour une paire de stations on peut avoir 9 paires de traces (stations différentes) ou 3 paires de traces (stations identiques) ?
    • prendre en compte les autocorrélations (même si on sait que le δt sera 0) : pour une paire de stations on peut avoir 9 paires de traces (stations différentes) ou 6 paires de traces (stations identiques) ?
    • ca biaise les proportions...
  • Figures films de carte de nombre de voisins résiduel en fonction du niveau de cohérence associée.
    • on confirme l'idée d'un seuillage sur le niveau de cohérence plutôt que sur l'interdistance
    • Films:
films de la carte (lon,lat) du niveau de cohérence pour chaque paire (un trait de couleur pour représenter une paire), pour itération 1 et 2, pour chacune des 9 composantes
films de la carte (lon,lat) du nombre de stations (ou de traces ?) qui restent connectées après application d'un seuil de cohérence donné

[edit] 19/09/2017 : Réunion de rentrée

Participants : G. Cougoulat, L. Métivier, C. Picard, P. Roux, P.-A. Bouttier, V. Louvet, L. Tavard, A. Lecointre

Durée : 2h

  • 1. Get inputs et accès NFS resif9@BIO depuis luke : point d'information sur le travail de Bruno et Pierre V.

il y avait plusieurs points à traiter:

- faire en sorte que les données soient accessibles à la fois depuis les machines ISterre et depuis les machines CIMENT, en tenant compte de la différence des référentiels utilisateur.

- tenir compte du fait que les usagers peuvent appartenir à plus d'un groupe projet (au sens POSIX), notamment dans CIMENT (=> les droits POSIX classiques ne sont pas suffisants pour exprimer cela).

- conserver la "privacy" des données.

Jusqu'à présent, l'approche avait été d'utiliser la fonction de remapping d'UID/GID dans NFS, mais cette technique n'a pas fonctionné de façon satisfaisante du côté CIMENT. Ce n'était de toute façon qu'une rustine, dont on aurait vite atteint les limites.

Bruno a suggéré de mettre en place des droits complémentaires sous formes d'ACL => c'est ce que nous avons fait, et supprimé le remapping d'UID. Cela a fonctionné côté CIMENT, mais par effet de bord nous avons perdu les accès côté ISterre. Suite à échange avec Rodolphe, le problème a été localisé dans la conf du serveur NFS sur resif9 (=> correction faite).

Donc : les modalités d'accès restent les mêmes depuis le cluster ISterre (groupe LDAP dédié) et les accès aux usagers CIMENT sont gérés par ajout d'ACL, au cas par cas. Cela reste donc souple.

  • 2. Get Input et Webservices

Suite réunion Glenn+Violaine+Albanne juillet → nécessité d’adapter le script OAR de « GetInputs+Preprocess » pour ne pas seulement attaquer les données de resif9@BIO en NFS mais être capable d’utiliser les webservices. Les 2 options supplémentaires :

- wget (GET method = 1 requête / daily trace)

- obspy.fdsn.client (1 exécution du code python / daily trace)

sont maintenant opérationnelles. wget semble un peu plus efficace.

C’était le prérequis pour étudier le workflow GetInputs+Preprocess sous l’approche Spark/Hadoop.

Note: la génération du fichier de paramètre stationByCodes.txt devrait aussi pouvoir se faire par requête webservices sur les métadata (réseau, lonmin lonmax, latmin, latmax) plutôt que via le portail web RESIF. Seulement là on aura portabilité totale d'un centre de données à un autre, et ca évitera les longues attentes de réparation du portail web (en carafe depuis le 31 août...). -> To Do -> Albanne

  • 3. Workflow prétraitements:

Approche Spark/Hadoop Système de fichiers distribués : point d'info discussion de cet été avec Violaine, Glenn, Christophe Cancé

- Compte Glenn au LAL.

- GRICAD possède quelques serveurs sur lesquels ils veulent installer des systèmes de fichiers distribués (SPARK/HADOOP) pour répondre aux besoins de formation d’une part et d’utilisation et d’expérimentation d’autre part. → réunion de lancement le 3 octobre.

Dans tout les cas : Portage du script OAR GetInputs+Preprocess à prévoir. -> session de travail prévue en semaine 25/09/2017 Glenn + Albanne

Empreinte mémoire du code de prétraitements : plusieurs centaines de MB (2GB si RAPATRIE_OPTION=0, 800MB si RAPATRIE_OPTION=3). Alors qu'une trace journalière 200Hz F64 avant décimation ne devrait pas dépasser 130MB. A explorer: P.-A. se charge d'une analyse mémoire et perfs du code à partir de la révision SVN 241.

  • 4. Prétraitements et la méthode de Laurent (enveloppe autour de plusieurs bandes de fréquence) : jamais testée encore...
  • 5. Calcul des δt

Les plots de δt [120,2,14535] et cohérences [120,2,14535] associées.

- Distribution des cohérences : fit gaussienne. Caractère aléatoire des évènements dans une loi normale. Indépendance ?

- Quelle est la signification des cohérences négatives (40 % de points auxquels on applique δt = ± T/2) ? :

Exemple figure FR.ESCA.01.BHZ_FR.TURF.00.BHZ 2017011 (m=60) H17 : on applique +T/2 suivi de -T/2 → pas d’effet au final ?

Exemple figure FR.ESCA.01.BHZ_FR.TURF.00.BHZ 2017012 (m=60) H09 : on applique -T/2 puis 0 → on a créé un δt qui n’existait pas au départ… Ici un seuil en cohérence n’empêcherait pas ce pb.

Dans les 2 cas m=120min n'a pas ce problème.

Action 1: Tester une autre méthode de calcul de δt pour ne plus être contraint par cet échantillonnage 0.1sec : faire le calcul du δt via le rapport des phases des deux signaux à une frequence choisie (la frequence centrale du spectre par exemple).

X=0.15; %frequence choisie
temps=0.1:0.1:...; % vecteur temps de meme longueur que vecteur signal [npc=6001]
phase_1=exp(2*pi*i*X*temps).';
DFT_base=mean(signal_base.*phase_1,1); % DFT signal ref
DFT_fluct=mean(signal__fluct.*phase_1,1); % DFT signal perturbé
deltat=angle(DFT_fluct./DFT_base)/2/pi/X; %calcul du dt

Action 2: Se restreindre aux paires de composantes verticales: donc 57*56/2 = 1596 paires (EXP6) (on garde pour plus tard cette idée d'améliorer le SNR en sommant sur les 9 composantes du tenseur car il faudrait préalablement corriger du déphasage les différentes composantes)

  • 6. Calcul des Δt.

- rappel de la méthode AtA=N*Id si hypothèse de fermeture pour obtenir les Δt à partir des δt → approche testée avec jeu de données Pyrope (avec succès) et avec jeu de données AlpArray (échec, on garde trop de déchets si on n'applique aucun seuil de cohérence).

- Idée: on n'applique pas l'hypothèse de fermeture.

Soit A[N,N(N-1)/2] matrice d'incidence du graphe.

Soit AT*A la matrice de Laplace (=représentant un graphe).

Soit B la matrice d'adjacence: ai,j; i .ne. j est le nombre d'arêtes liant le sommet i au sommet j. ai,i est le nombre de boucle au sommet i (toujours nul pour nous).

Soit D la matrice diagonale des degrés de chaque sommet du graphe (nombre de liens aboutissant à ce sommet).

Note: AT*A, B, (et bien sur D) sont toujours symétriques car notre graphe est orienté.

AT*A = D - B , donc ce serait intéressant de pouvoir diagonaliser B ?

...

Travail en cours, -> Christophe

- normalisation finale par rapport à une trace dont on sait que Δt=0 (car jamais perdu GPS) : cf Sens-Schönfelder : différent de notre hypothèse de fermeture. Quelle station pour ce rôle de référence absolue ? à rediscuter plus tard.

  • 7. Une présentation aux journées techniques et scientifiques RESIF ?
  • 8. Meeting hors les murs à Chamonix : Dates fixées aux LU-MA 18-19 déc.
  • 9. ~12kE pour équipement (noeud luke ?, bettik ?, ...)
  • 10. Prospectives: Déploiement réseau denses Resolve (glacier d'Argentière, 100 capteurs très denses) -> Application de la méthode pour estimer les dérives spatiales à partir des Δt.

[edit] 19/09/2017 : Réunion informelle interne à ISTerre avec les responsables techniques du projet AlpArray et les responsables techniques RESIF

Participants : A. Paul, C. Aubert, D. Wolyniec, C. Pequegnat, A. Lecointre

Durée : 2h

L'objectif de cette présentation est de présenter et discuter du projet iWORMS avec les responsables techniques et scientifiques du projet AlpArray (C. Aubert, A. Paul) et RESIF (C. Pequegnat, D. Wolyniec).

  • Il faudrait explorer la relation entre distribution des cohérences et orientation azimutale des paires de stations considérées. En hiver (période d'étude jan 2015), la source de bruit sismique ambiant est localisé dans les GIN Seas (trains d'ondes de tempêtes hivernales à l'approche du relief bathymétrique). Les paires qui reconstruisent le mieux les corrélations (et pour lesquelles on peut espérer avoir de bonnes cohérences pour le calcul du δt) devraient être alignées dans la direction du bruit sismique ambiant (donc direction NW - SE pour notre subset AlpArray du quart SE de la France).
coherenceiter1 = f (azimuth)
2017012
iter1
Az.png


  • Une dépendance avec la saison est également attendue (en été la source de bruit sismique ambiant vient plutôt de l'Océan Indien).
  • On pourrait stacker davantage pour le calcul de notre GF de référence (exemple les 6 derniers mois glissants plutot que 24h). Mais pas dans le cas d'une dérive δt avérée (-> "serpent se mord la queue")
  • A tester: David fournit une liste de stations / dates ou une dérive Δt a effectivement été observée : il s'agit de station Z3 ayant "perdu" leur GPS:
    • A159A de 2017-06-11T00:00:00 à 2017-06-23T13:00:00 et de 2017-07-13T14:40:00 à 2017-07-22T11:00:00
    • A162A de 2017-07-26T07:00:00 à 2017-08-09T10:00:00
    • A198A de 2017-08-03T10:20:00 à 2017-08-09T12:30:00
A159A, A162A, A198A
A159A A162A A198A.png
  • A tester: Coralie fournit une station avec une dérive connue:
    • OGRV : en attente de la période
  • A tester: David et Catherine fournissent l'indication d'une station avec une mauvaise polarité:
    • FR.ESCA.01.HNE entre le 2014-07-30 et le 2017-06-26 08h50 - 08h55.
ESCA, MON, CALF
ESCA MON CALF.png
FR.ESCA.01.HNE-FR.MON.00.HNE
2017001 (2017-01-01) 2017174 (2017-06-23) 2017180 (2017-06-29) 2017244 (2017-09-01)
it1:
Cc iter1 FR.ESCA.01.HNE FR.MON.00.HNE 2017001.png
Cc iter1 FR.ESCA.01.HNE FR.MON.00.HNE 2017174.png
Cc iter1 FR.ESCA.01.HNE FR.MON.00.HNE 2017180.png
Cc iter1 FR.ESCA.01.HNE FR.MON.00.HNE 2017244.png
it2:
Cc iter2 FR.ESCA.01.HNE FR.MON.00.HNE 2017001.png
Cc iter2 FR.ESCA.01.HNE FR.MON.00.HNE 2017174.png
Cc iter2 FR.ESCA.01.HNE FR.MON.00.HNE 2017180.png
Cc iter2 FR.ESCA.01.HNE FR.MON.00.HNE 2017244.png
FR.ESCA.01.HNE - other traces
Coh FR.ESCA.01.HNE.png
FR.ESCA.01.HNE-FR.CALF.00.HNE
2017175 (2017-06-24) 2017178 (2017-06-27)
it1:
Cc iter1 FR.CALF.00.HNE FR.ESCA.01.HNE 2017175.png
Cc iter1 FR.CALF.00.HNE FR.ESCA.01.HNE 2017178.png
it2:
Cc iter2 FR.CALF.00.HNE FR.ESCA.01.HNE 2017175.png
Cc iter2 FR.CALF.00.HNE FR.ESCA.01.HNE 2017178.png



  • Etape finale de recalage Δt sur un temps absolu (G SSB serait une bonne candidate)
  • Journées techniques et scientifiques RESIF 10-12 Octobre Vendée : intervention en atelier challenges/enjeux : format (MSeed3/HDF5), workflows.

[edit] 25/09/2017 : Réunion Δt Pyrope

Participants : L. Stehly, P. Roux, A. Lecointre

Durée : 1h15

[Voir la section dédiée]

On va se focaliser sur la bande de fréquence 3-5sec.

On veut comparer les stack 1day aux stack 5days : impact de l'application du filtre avant (= sur les xcorr) vs après (= sur les Δt) le calcul des δt et Δt.

Les 6 stations qui présentent des Δt brutaux (PF01, PF06, PF03, PF09, PF12, PF15) sont des Tellus, ce matériel est connu pour avoir des dérives brutales d'1sec corrélées avec les journées de visite sur le terrain.

On veut rajouter des plots 2D de Δt = f(t), un plot pour chaque station (plus lisible que les imagesc).

On veut comparer aux plots de Laurent avec la méthode de la médiane (superposer les plots 2D).

Note: pour la calcul des δt, Laurent a fenêtré les xcorr autour de l'ondes de surface + 2 ou 3 périodes puis a rééchantillonné à 10x, 100x (?). Donc une méthode différente de celle de Sens-Schönfelder que l'on utilise en 1ère approche pour AlpArray.

Quelle info à postériori sur la fiabilité (sur la cohérence) ? Voir dt.coefficient dans les fichiers de Laurent. On remarque 30% de coefficients nuls associés à des δt nuls : seuillage ? Donc le calcul des δt aurait lui aussi été délicat pour PE/PF (30% de cohérence inférieures au seuil ?).

On veut comparer l'approche relative (méthode de la médiane pour δt -> Δt) et l'approche complète (ATA=N*Id avec Hypothèse de fermeture et nécessité d'un graphe complet) -> superposition de résultats Δt avec errorbar associées.

Filtre de Wiener de Ludovic Moreau : ajout au preprocess quand le calcul du δt s'avère délicat ?

Et si on forçait Σ δt = 0 avant calcul des Δt ?

[edit] 06/11/2017 : Réunion 10h 258 IMAG

Participants : L. Stehly, P. Roux, P.-A. Bouttier, B. Bzeznik, C. Picard, L. Métivier, A. Lecointre

Durée : 1h45

Compte-rendu:

  • 1. Workflow prétraitements:
    • Approche MapReduce Hadoop/Spark :
      • Idée: un framework intégrant : le //isme, la tolérance aux pannes (checkpointing, réplication), le déploiement
      • Idée: un outil qui continue à scaler (même si ce n'est pas la version la plus optimisée) qqsoit la quantité de données à traiter...
      • Actuellement on ne dispose pas de machines Hadoop (en cours d'installation mais difficultés OS / carte réseau)
      • Glenn et Albanne suivent la formation M2-GI-ENSIMAG (V. Leroy) Hadoop/MapReduce (TP en Java) + Spark/Stream (TP en Scala) en nov-déc 2017 (sur machines ENSIMAG dédiées aux formations)
      • les "keys" du //isme "Map" seraient les périodes: on aimerait tester l'effect de "keys" horaires/journalières/mensuelles sur la scalabilité. (actuellement avec le //isme "grille" la granularité est sur les jours sans option de modification).
    • Preprocess code: analyse empreinte mémoire. Le code doit être modularisé pour être profilé (Pierre Antoine s'en charge ainsi que du profilage). Pierre-Antoine et Albanne se verront le 13 novembre pour :
      • précisions sur le code
      • explications sur les inputs pour s'assurer une exécution rapide pour les tests (=reduire la boite lat/lon)
      • Albanne doit auparavant modifier le script de soumission pour remplacer l'input stationsByCodes.txt (issu du portail web resif) par une requête webservices sur les metadata (requête "Station") pour faciliter cette gestion des paramètres lat/lon (cette évolution était de toute façon prévue et motivée par Catherine car le portail web RESIF est trop souvent indisponible et nos inputs ne peuvent pas se baser dessus).
  • 2. Calcul des δt:
    • [quel lag pour les corr : 300sec -> reduire à 150sec ?] Effet: A la 1ère itération, le pic de coherence positive se déplace de 0.30-0.35 vers 0.40 (idem pour les cohérence négatives). Donc réduire nn_lag ne résoudra pas le problème. Il faut fenêtrer sur l'onde de surface à partir d'un stack sur "une grande période" et manipuler cette info (limite_fenetre=f(paire de traces)) qui doit devenir un input du code -> Laurent fournit un script MatLab qui décrit la méthode pour fenêtrer sur l'onde de surface.
    • Pour répondre au problème de l'échantillonnage 0.1sec : on propose de [faire le calcul du δt via le rapport des phases des deux signaux à une frequence choisie (la frequence centrale du spectre par exemple)]. Tant qu'on n'a pas fenêtré sur l'onde de surface, ca ne marche pas mieux. Il faut qu'on fournisse deux versions de code : une avec nn_lag ajustable (on l'a déjà), une avec fenétrage sur l'onde de surface (code à écrire et caractéristiques des fenêtres pour chaque paires de traces à créer). Ca pose un certain nombre de questions:
      • il faut fabriquer cet input à partir de données sur une longue période (1an serait souhaitable), toutes nos paires ne sont pas forcément valides sur cette même période d'un an
      • on part d'un input (=une hypothèse de départ) basé sur un an de données qui sont susceptibles d'être affectées par une erreur de temps
      • on doit définir une fenêtre pour chaque paire de trace, et pas simplement pour chaque paire de station (?)
      • quand ca ne marche pas même en stackant sur un an (onde de surface indétectable), que fait on ?
      • ...
  • 3. δt -> Δt:
    • File:Dt Dt.pdf : point avec Christophe et tous : travail en cours ...
    • on confirme que Σ Δt ≠ 0
    • on a Σ δt = 0 le long du plus grand circuit, c'est cette hypothèse qu'il faut utiliser ainsi que la réorientation de certains circuits pour réécrire G (?)
    • Dans l'immédiat, on utilisera la méthode de la médiane pour évaluer les expériences prévues au point 2 ci dessus.
  • 4. Meeting iWORMS Chamonix 18-19déc ([8pers/12])
    • Téléphérique aiguille du midi benne à 11h donc départ campus vers 8h
    • programme scientifique du mardi ?
    • horaires départ/retour ?
    • picnic lundi midi ?
    • voitures ?

[edit] 29/11/2017 : Réunion 10h 206 IMAG

Participants : Pierre-Antoine Bouttier, Bruno Bzeznik, Catherine Pequegnat, Violaine Louvet, Christophe Picard, Philippe Roux, Laurent Stehly, Albanne Lecointre

Durée : 2h

Ordre du jour:

  • code Python prétraitements
  • tests de montée en charge wget sur RESIF@SUMMER ?
  • δt -> Δt
  • Proposition de programme de la demie journée de travail du 19/12/2017 à Chamonix:
    • 9h:
    • La corrélation de bruit et les applications en géophysique
    • iworms revue de projet
    • δt -> Δt
    • quels workflows pour quelles infrastructures ?
    • 12h30 : repas puis retour

[edit] 18-19/12/2017 : Réunions Chamonix

Participants : Pierre-Antoine Bouttier, Bruno Bzeznik, Catherine Pequegnat, Violaine Louvet, Philippe Roux, Christophe Picard, Albanne Lecointre

Programme et Compte-Rendu :

  • 18/12/2017 14:30-15:30 : Contexte scientifique (animateur P. Roux) File:Chamonix.pdf
  • 18/12/2017 17:00-18:00 : Le projet HPCDA@UGA (animateur B. Bzeznik) File:Dahu.pdf
  • 19/12/2017 9:00-10:00 : iWORMS méthodologies de calcul des erreurs d'horloge (animateur C. Picard) File:Gdt.pdf

[edit] 23/01/2018 : Réunion 10h 248 IMAG

Participants : Pierre-Antoine Bouttier, Violaine Louvet, Philippe Roux, Christophe Picard, Glenn Cougoulat, Albanne Lecointre

Durée : 1h30

Compte-Rendu:

  • Prétraitements:
    • Mémoire
    • //istaion multicoeurs joblib (et multinoeuds dask ?)
    • Tests montée en charge WGET SUMMER
  • δt et Δt: Mise en place de tests avec données synthétiques
  • papier Sens-Schönfelder: quelles différences méthodologiques par rapport à nous : cf slides 10 à 14 dans mes slides Chamonix
  • figures distribution = f(dt,coherence) sans la normalisation:
2017012 - timestep=1h - iteration1 - zoom sur onde de surface (côté enveloppe max)
dt mesuré par la corrélation entre GF1h et GF1d, et cohérence associée dt mesuré par le rapport des phases entre GF1h et GF1d, et niveau de cohérence associée à l'amplitude de la corrélation pour ce dt (le plus proche à 0.1s près)
FR.OGGM.00.HHZ-FR.OGSM.00.HHZ
Cc iter1 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 sens zoomSWcotemax.png
Cc iter1 FR.OGGM.00.HHZ FR.OGSM.00.HHZ 300 phase zoomSWcotemax.png
2346 HHZ pairs
Dt distribution 2017012 iter1.absSW P nonorm.png
Dt distribution 2017012 iter1.absSW P phase coh nonorm.png

[edit] 26/03/2018 : Réunion 15h 248 IMAG : Mesure des δt sur des GF synthétiques complètes bruitées mais synchronisées; Calcul des Δt associées

Participants : P. Roux, C. Picard, A. Lecointre

Durée : 2h15

Compte-rendu :

0) Technique de détection de phase dans le domaine réel Ecriture d'un filtre spectral pour détecter ces changements de phase ? Oublié: à reprogrammer pour la prochaine réunion

1) Prétraitements: Annulé (Pierre-Antoine absent): [réunion de travail programmée 28/03] (P.-A.B. + A.L.) pour (re)définir plan d'actions

2) Mesure des δt: Mise en place de tests avec données synthétiques: création de jeux de données synthétiques (GF) et mesure de δt associés. (Voir code ist-oar:/home/lecoinal/IWORMS_tools/synthetic_GF/synthetic_GF.m)

  • jeu de données synthétiques dans lesquels on n'a simulé aucun décalage, toutes les stations sont synchronisées
  • une seule composante HHZ
  • 4 réseaux: FR, Z3, G, ZH
  • box: lat 43N - 46N ; lon 3E - 8E
  • 1 year 2016183 -> 2017181
  • -> 79 traces -> 79*78/2=3081
lecointre@luke4:/scratch_md1200/lecointre/IWORMS/METRICS/SW$ wc -l pair_dist.txt stalatlon.txt
 3081 pair_dist.txt
   79 stalatlon.txt

Méthode pour créer des GF synthétiques:

  • double gauspuls 10Hz, symétrique par rapport à 0, sur lags -300sec:300sec [6001]. la distance du puls au zero est donnée par la distance entre les 2 stations.
  • une matrice r[24,6001] de random noise sur range [-1 ... 1]
  • filtre chaque ligne de r entre 0.1 et 0.2Hz
  • r=r/max(abs(r));  % normalize
  • r=r/SNR;
  • ajoute chacun des 24 r au double gauspuls : ce sont nos 24 GF synthétiques pour une paire
  • moyenne des 24: ce sera la GF de référence pour cette paire
  • répète l'opération pour les 3081 paires

Il nous faut les distances pour chacune des 3081 paires (pair_dist.txt) mais aussi la position de l'onde de surface supposée dans la GF[-300sec:+300sec]. On estime cette position (code SW_guess_from_dist.m):

t1 = dist / 5                   ! [dist en km et t en sec]
t2 = dist / 1.5
dt = 0.1                        ! 10Hz
if t2-t1 < dt then t2 = t1+5*dt ! can happend for close stations
if t1 > 300 then Abort          ! distance trop grande, en dehors de la fenetre de correlation (ca n'arrive pas avec nos distances)
t1 = MAX(t1,0)                  ! cut the [t1,t2] inside [0,300] window (ca n'arrive pas par construction)
t2 = MIN(t2,300)                ! (ca n'arrive pas avec nos distances)

On obtient une liste:

$ head -5 dist_SW_synthetiques.dat
39760.246 3081 3266 % FR.ARBF.00.HHZ FR.ARTF.00.HHZ
135732.047 3272 3906 % FR.ARBF.00.HHZ FR.BANN.00.HHZ
76790.758 3155 3513 % FR.ARBF.00.HHZ FR.BLAF.00.HHZ
42586.914 3086 3285 % FR.ARBF.00.HHZ FR.BSTF.00.HHZ
131300.703 3264 3876 % FR.ARBF.00.HHZ FR.CALF.00.HHZ
$ sort -nk 1 dist_SW_synthetiques.dat|head -1
0.000 3001 3006 % FR.OGS1.00.HHZ FR.OGS2.00.HHZ
$ sort -rnk 1 dist_SW_synthetiques.dat|head -1
409177.438 3819 5729 % FR.CFF.00.HHZ FR.MON.00.HHZ

Remarques:

  • Toujours du côté positif (ca n'a de toute façon pas d'importance pour ces gauspuls synthétiques qui sont symétriques)
  • La largeur de la fenêtre varie avec la distance. Cette largeur est comprise entre 5*dt=0.5sec (seulement 6 points !) pour les paires les plus proches et 191sec (1911 points) pour la paire la plus éloignée (409km). Ca révèle un biais potentiel sur les niveaux de cohérence qu'il ne faudra pas interpréter comme un effet de la distance. On pourrait envisager une autre méthode pour déterminer la largeur de fenêtre: peut être une largeur qui augmente avec la distance comme c'est le cas actuellement (dans ce cas est on d'accord avec la largeur min = 0.5sec = 6 échantillons seulement ?), mais avec un certain seuil (quelle largeur seuil ?, 40sec = 401 échantillons comme le code cc2da.py ?).
  • Les limites de vitesses ont été fixées (arbitrairement) à 1.5km/s et 5km/s. On pourrait envisager une autre méthode: hypothèse milieu isotrope, une section sismique (de nos GF synthétiques ou moyennée sur un an si données réelles), mesure de la pente pour extraire une vitesse moyenne (dans le cas synthétique on retrouverait logiquement vel=3km/s qui est le paramètre utilisé pour le décalage de nos gauspuls synthétiques), détermination d'un fenêtre de largeur fixe autour de cette vitesse moyenne.
  • On aurait pu réutiliser les positions SW qui nous avions estimé (voir [fenêtrage sur l'onde de surface]) à partir des données réelles (moyennes annuelles des GF) pour au moins les 2838 paires de données réelles (sur un an) parmi ces 3081 paires synthétiques.
    On ne le fait pas car certaines paires de données réelles, même en moyennant sur un an, ne permettent pas d'obtenir une GF qui permette d'estimer correctement sa position.
    Exemple:
128500.688 3580 3980 % FR.RSL.00.HHZ Z3.A186B.00.HHZ
analyse avec position SW déduite de données réelles position SW déduite de données réelles analyse avec position SW synthétique
SNR 100 SW 1 dist 128.5007.png
Fig 1969.png
SNR 100 SW 1 dist 128.5007 distSWsynt.png
Ici si on utilise la position déduite de la SW réelle, on vise à côté du gauspuls synthétique et on obtient des δt ≠ 0 associés à des cohérences très faibles En pointillés la position SW estimée à partir de la distance uniquement (donc qui serait centrée sur gauspuls synthétique). En traits plein la position SW retenue si on ajoute l'information de la SW réelle moyenne sur un an (qui apparaît mal reconstruite ici). Ici si on utilise la position déduite uniquement de la distance, on vise bien sur le gauspuls et on obtient bien δt = 0 et coh = 1.

Exemples détaillés de mesure de δt pour ipair=1651 (distance=0km), ipair=1 (distance 40km) et pour ipair=525 (374km):

On présente ci dessous des exemples pour lesquels la méthode fonctionne bien, que l'on effectue la mesure sur la fenêtre complète [-300:300sec] ou sur la fenêtre zoomée sur SW (largeur de fenêtre 0.5sec ou 35sec ou 170sec selon la distance)

On effectue la mesure du δt avec la méthode rapport des phases.

Sans surprise, (dans les cas SNR=100 ou SNR=10), les niveaux de cohérence sont supérieurs (et les δt sont encore plus faibles) quand on zoome sur l'onde de surface.

On a l'impression que c'est le cas SNR=1 qui s'apparente le plus à nos données réelles (voire aux meilleures de nos données réelles).

lecoinal@ist-oar:~/IWORMS_tools/synthetic_GF$ head -1 dist_SW.dat
39760.246 3001 3336 % FR.ARBF.00.HHZ FR.ARTF.00.HHZ
lecoinal@ist-oar:~/IWORMS_tools/synthetic_GF$ grep -n "FR.CFF.00.HHZ ZH.POSAS.00.HHZ" dist_SW.dat
525:374127.719 3749 5495 % FR.CFF.00.HHZ ZH.POSAS.00.HHZ
lecoinal@ist-oar:~/IWORMS_tools/synthetic_GF$ sort -nk 1 dist_SW_synthetiques.dat|head -1
0.000 3001 3006 % FR.OGS1.00.HHZ FR.OGS2.00.HHZ
lecoinal@ist-oar:~/IWORMS_tools/synthetic_GF$ grep -n "FR.OGS1.00.HHZ FR.OGS2.00.HHZ" dist_SW_synthetiques.dat
1651:0.000 3001 3006 % FR.OGS1.00.HHZ FR.OGS2.00.HHZ


(random noise bandpass 5-10sec avec filtfilt)
0km, [-300:300sec] (filter) 0km, zoom sur SW (filter et SWwindow synthétique) 0km, zoom sur SW fenetre largeur fixe 40sec 0km, zoom sur SW fenetre largeur fixe 40sec et decrease GF with distance 40km, [-300:300sec] 40km, zoom sur SW 40km, zoom sur SW fenetre largeur fixe 40sec 40km, zoom sur SW fenetre largeur fixe 40sec et decrease GF with distance 374km, [-300:300sec] 374km, zoom sur SW 374km, zoom sur SW fenetre largeur fixe 40sec 374km, zoom sur SW fenetre largeur fixe 40sec et decrease GF with distance
SNR=100
SNR 100 SW 0 dist 0.png
SNR 100 SW 1 dist 0.png
SNR 100 SW 2 dist 0.png
SNR 100 SW 2 dist 0 decGF.png
SNR 100 SW 0 dist 39.7602.png
SNR 100 SW 1 dist 39.7602.png
SNR 100 SW 2 dist 39.7602.png
SNR 100 SW 2 dist 39.7602 decGF.png
SNR 100 SW 0 dist 374.1277.png
SNR 100 SW 1 dist 374.1277.png
SNR 100 SW 2 dist 374.1277.png
SNR 100 SW 2 dist 374.1277 decGF.png
SNR=10
SNR 10 SW 0 dist 0.png
SNR 10 SW 1 dist 0.png
SNR 10 SW 2 dist 0.png
SNR 10 SW 2 dist 0 decGF.png
SNR 10 SW 0 dist 39.7602.png
SNR 10 SW 1 dist 39.7602.png
SNR 10 SW 2 dist 39.7602.png
SNR 10 SW 2 dist 39.7602 decGF.png
SNR 10 SW 0 dist 374.1277.png
SNR 10 SW 1 dist 374.1277.png
SNR 10 SW 2 dist 374.1277.png
SNR 10 SW 2 dist 374.1277 decGF.png
SNR=1
SNR 1 SW 0 dist 0.png
SNR 1 SW 1 dist 0.png
SNR 1 SW 2 dist 0.png
SNR 1 SW 2 dist 0 decGF.png
SNR 1 SW 0 dist 39.7602.png
SNR 1 SW 1 dist 39.7602.png
SNR 1 SW 2 dist 39.7602.png
SNR 1 SW 2 dist 39.7602 decGF.png
SNR 1 SW 0 dist 374.1277.png
SNR 1 SW 1 dist 374.1277.png
SNR 1 SW 2 dist 374.1277.png
SNR 1 SW 2 dist 374.1277 decGF.png


Distribution des δt mesurés

On s'intéresse à la distribution des ABS(δt) en fonction du SNR et du zoom SW dans les 3 cas suivants:

  • mesure effectuée sur la fenêtre [-300:300]sec tout entière (colonne 1)
  • mesure effectuée sur la fenêtre de l'onde de surface telle que déterminée à partir des infos (distance,GFannuelle de données réelles) (colonne 2 : points rouge et noir)
  • mesure effectuée sur la fenêtre de l'onde de surface telle que déterminée à partir des infos (distance uniquement) (colonne 2 : points cyan)

On compare également le type de filtre bandpass 0.1-0.2Hz appliqué au random noise:

  • filter MatLab ordre 3 (points noirs et cyan)
  • filtfilt MatLab ordre 1 (zero phase filtering) (points rouge)
distribution δt
[-300:300sec] zoom sur SW
SNR=100
Dt dist SNR 100 SW0.png
Dt dist SNR 100 SW1.png
SNR=10
Dt dist SNR 10 SW0.png
Dt dist SNR 10 SW1.png
SNR=1
Dt dist SNR 1 SW0.png
Dt dist SNR 1 SW1.png

Les courbes rouges (filtfilt) sont plus piquées que les courbes filter. A rediscuter avec Philippe...

On normalise les plots ci dessus (normalise par trapezoidal integration) et on vérifie que l'on fitte bien une gaussienne (distinguer + et -), puisqu'on a ajouté du bruit gaussien:

distribution δt (mesure δt avec méthode rapport des phases)
[-300:300sec] zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) et exclude points dt=0 avant fit zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec (401points) centrée sur dist/V=dist/3km/s (sauf si à cheval sur + et -) zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec (401points) centrée sur dist/V=dist/3km/s (sauf si à cheval sur + et -) et décroissance de GF amplaitude avec la distance
SNR=100
Dt dist gaussianfit SNR 100 SW0.png
Dt dist gaussianfit SNR 100 SW1.png
Dt dist gaussianfit SNR 100 SW1 exclude.png
Dt dist gaussianfit SNR 100 SW2.png
Dt dist gaussianfit SNR 100 SW2 decGF1.png
SNR=10
Dt dist gaussianfit SNR 10 SW0.png
Dt dist gaussianfit SNR 10 SW1.png
Dt dist gaussianfit SNR 10 SW1 exclude.png
Dt dist gaussianfit SNR 10 SW2.png
Dt dist gaussianfit SNR 10 SW2 decGF1.png
SNR=1
Dt dist gaussianfit SNR 1 SW0.png
Dt dist gaussianfit SNR 1 SW1.png
Dt dist gaussianfit SNR 1 SW1 exclude.png
Dt dist gaussianfit SNR 1 SW2.png
Dt dist gaussianfit SNR 1 SW2 decGF1.png


Distribution des δt mesurés avec info du niveau de cohérence associée

On s'intéresse à la distribution des δt et coherence associée en fonction du SNR et du zoom SW ou non:

Ici on ne plotte que le cas du bandpass (methode filter et non filtfilt) car on ne voit pas d'impact sur ces plots.

distribution δt et coherence
[-300:300sec] + mesure δt avec méthode rapport des phases [-300:300sec] + mesure δt avec méthode corrélations(GF_ref,GF_timestep) zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) + mesure δt avec méthode rapport des phases zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) + mesure δt avec méthode corrélations(GF_ref,GF_timestep) zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec ou 401pts (uniquement = f(dist,V=3km/s)) + mesure δt avec méthode rapport des phases zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec ou 401pts (uniquement = f(dist,V=3km/s)) + mesure δt avec méthode corrélations(GF_ref,GF_timestep) zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec ou 401pts (uniquement = f(dist,V=3km/s)) + mesure δt avec méthode rapport des phases + decGF=1
SNR=100
Dt distribution SNR 100 SW0 FILTER.png
Dt distribution SNR 100 SW0 CC1 CC FILTER.png
Dt distribution SNR 100 SW1 FILTER SW SYNT.png
Dt distribution SNR 100 SW1 CC1 CC FILTER.png
Dt distribution SNR 100 SW2 CC0 FILTER.png
Dt distribution SNR 100 SW2 CC1 CC FILTER.png
Dt distribution SNR 100 SW2 CC0 FILTER decGF.png
SNR=10
Dt distribution SNR 10 SW0 FILTER.png
Dt distribution SNR 10 SW0 CC1 CC FILTER.png
Dt distribution SNR 10 SW1 FILTER SW SYNT.png
Dt distribution SNR 10 SW1 CC1 CC FILTER.png
Dt distribution SNR 10 SW2 CC0 FILTER.png
Dt distribution SNR 10 SW2 CC1 CC FILTER.png
Dt distribution SNR 10 SW2 CC0 FILTER decGF.png
SNR=1
Dt distribution SNR 1 SW0 FILTER.png
Dt distribution SNR 1 SW0 CC1 CC FILTER.png
Dt distribution SNR 1 SW1 FILTER SW SYNT.png
Dt distribution SNR 1 SW1 CC1 CC FILTER.png
Dt distribution SNR 1 SW2 CC0 FILTER.png
Dt distribution SNR 1 SW2 CC1 CC FILTER.png
Dt distribution SNR 1 SW2 CC0 FILTER decGF.png

Remarque : Quand SNR=100 et SW=0 ([-300:300sec] + méthode rapport des phases), on observe des points avec faible cohérence et δt aberrant. Pourquoi ?

On a vérifié que ce n'est pas le filtre filtfilt vs filter.

On dirait un effet de la méthode rapport des phases par rapport à la méthode corrélation. Voir script jouet example_badcoh.m, ou l'on reproduit facilement , sur [-300:300sec], des exemples ou :

  • méthode rapport des phases --> δt élevé avec cohérence faible
  • méthode corrélation --> δt=0 avec coh=1

Exemple:

Dt coh distance SNR100 SW0 CC0.png
Phase vs cc SNR100 SW0 ipair2279.png
Phase vs cc SNR100 SW0 ipair1.png
  • C'est toujours pour certaines distances que le problème se pose (par exemple 134km ci dessus)
  • A 134km, les spectres sont bien identiques, mais s'annulent tout deux à la fréquence centrale 0.15Hz


Evolution du δt et de la cohérence associée au cours des 24 pas de temps. Les indices des paires en abscisse sont triées en distance croissante.

Evolution δt et cohérence (filter), mesure δt avec méthode rapport des phases
zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec (f(dist,V=3km/s)) zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec (f(dist,V=3km/s)) et decGF=1
SNR=100
Dt SNR 100 SW1 FILTER SW SYNT.png
Dt SNR 100 SW2.png
Dt SNR 100 SW2 decGF1.png
SNR=10
Dt SNR 10 SW1 FILTER SW SYNT.png
Dt SNR 10 SW2.png
Dt SNR 10 SW2 decGF1.png
SNR=1
Dt SNR 1 SW1 FILTER SW SYNT.png
Dt SNR 1 SW2.png
Dt SNR 1 SW2 decGF1.png

On remarque une dépendance à la distance. Il s'agit d'un biais lié aux largeurs variables des fenêtres SW. Attention à ne pas faire de fausse interprétation quand on travaillera avec des données réelles.

Les δt ne sont pas indiqués en sec mais en nombre de période (δt / Tmoyen).

  • SNR=100 : On mesure des δt de l'ordre de 0.03/7.5=0.004 périodes. Bref c'est 0.
  • SNR=10  : 4% de période
  • SNR=1  : 40% de période
  • TODO: ajouter les plots du cas de fenêtres de largeur fixe au delà d'un certain seuil

Autres idées en attente:

  • On a fait une référence 24 timesteps. Test aussi la référence 365 timesteps ?
  • Info distance pour de futurs jeux de données synthétiques ? en 1/sqrt(r) ...?
  • Prévoir futurs jeux de données synthétiques dans lesquels on simulera un décalage de la station A par rapport à toutes les autres ou un décalage de 2 stations par rapport à toutes les autres.

3) Calcul des Δt à partir de ces jeux de données de δt synthétiques.

On a donc construit des (12) jeux de données (BPfiltre={filter,filtfilt);window={SurfaceWaveZoom,[-300:300sec]};SNR={100,10,1}) de δt "synthétiques".

Exemple:

lecoinal@ist-oar:~/IWORMS_tools/synthetic_GF/FILTFILT$ h5ls -v SNR_100_SW1_rapportdesphases.h5
Opened "SNR_100_SW1_rapportdesphases.h5" with sec2 driver.
coh                      Dataset {3081/3081, 24/24}
   Location:  1:595568
   Links:     1
   Chunks:    {3081, 24} 591552 bytes
   Storage:   591552 logical bytes, 591552 allocated bytes, 100.00% utilization
   Type:      native double
dt_bypair                Dataset {3081/3081, 24/24}
   Location:  1:800
   Links:     1
   Chunks:    {3081, 24} 591552 bytes
   Storage:   591552 logical bytes, 591552 allocated bytes, 100.00% utilization
   Type:      native double

Note: datasets à transposer pour calcul Δt

Ces jeux de données synthétiques sont sur un graphe complet (nouveauté par rapport aux données réelles). On pourrait appliquer la méthode décrite dans les slides de Christophe de décembre 2017.

Revoir slide 8 ensemble:

  • c'est bien m(1)=cste (hypothèse sur le modèle) et non pas b(1)=cste (pas sur les observations).
  • on clarifie les notations: ne pas confondre d (les observations) vs d (les termes diagonaux de D la matrice diagonale issue de la diagonalisation du laplacien du graphe).
  • sera t il possible de prendre plutot m(k) avec k=2,N mais pas k=1 (pour prévoir le cas ou on aura imposé m(1)=Δt(1) = un décalage par rapport à toutes les autres stations.

Exemple des données synthétiques N traces (N=79):

  • On supose Δt(N)=cste=0
  • La solution devient unique, on peut la déterminer analytiquement avec un pivot de Gauss, elle est : pour k=1,...,N-1: Δt(k) = 1/N ( ∑ δtk,i; i=1,N; i≠k + ∑ δti,N; i=1,N-1)
    Algoritmiquement, cela revient à réarranger le vecteur de δt[N(N-1)/2] selon une matrice[N,N] antisymétrique avec des 0 sur la diagonale. Puis à effectuer une somme selon les lignes (ou les colonnes peu importe). Notre dataset de dimension δt[N(N-1)/2] se réduit donc à un dataset de dimension S[N]. Chaque Δti s'obtient par la différence (1/N)*(Si-SN) (si N est choisie comme référence).
  • TODO: Vérifier que l'on retrouve cette formule en passant par la décomposition en valeurs singulières, pour le cas N=4 par exemple.
  • Pour les cas SW={0,1} et SNR={100,10,1}, on obtient:
Si=∑jδti,j et Δti si l'on suppose que N=79 est référence (ΔtN=79=0) (méthode filter pour random noise ; mesure δt avec méthode rapport des phases)
Si ([-300:300sec]) Δti ([-300:300sec]) Si(zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s))) Δti (zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s))) Si(zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec) Δti (zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec) Si(zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec) et decGF=1 Δti (zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec) et decGF=1
SNR=100
Dt SNR100 SW0 sumdt.png
Dt SNR100 SW0.png
Dt SNR100 SW1 sumdt.png
Dt SNR100 SW1.png
Dt SNR100 SW2 sumdt.png
Dt SNR100 SW2.png
Dt SNR100 SW2 decGF1 sumdt.png
Dt SNR100 SW2 decGF1 ref79.png
SNR=10
Dt SNR10 SW0 sumdt.png
Dt SNR10 SW0.png
Dt SNR10 SW1 sumdt.png
Dt SNR10 SW1.png
Dt SNR10 SW2 sumdt.png
Dt SNR10 SW2.png
Dt SNR10 SW2 decGF1 sumdt.png
Dt SNR10 SW2 decGF1 ref79.png
SNR=1
Dt SNR1 SW0 sumdt.png
Dt SNR1 SW0.png
Dt SNR1 SW1 sumdt.png
Dt SNR1 SW1.png
Dt SNR1 SW2 sumdt.png
Dt SNR1 SW2.png
Dt SNR1 SW2 decGF1 sumdt.png
Dt SNR1 SW2 decGF1 ref79.png
  • C'est étrange qu'on ait des lignes toutes bleues ou toutes rouges. La probabilité d'en avoir une est d'1/79: en effet on avait une probabilité p=1/79 de choisir comme station de référence celle qui va être systématiquement en avance ou en retard par rapport à toutes les autres. C'est étrange d'en avoir plusieurs pour seulement 24 pas de temps. Il faut vérifier si pour ces pas de temps, les δt79,i étaient bien tous du même signe ? NON en fait si tout les Δti sont du même signe, cela signifie juste que SN était un max ou un min parmi tous les Si, ce que l'on vérifie sur les figures ci dessus (les symboles noirs représentent les min,max,median de Si pour chaque pas de temps)
    • je me demande s'il ne faudrait pas choisir comme référence i tel que Si soit mediane ? mais on ne veut pas changer de référence à chaque pas de temps ?
  • A priori on peut changer la référence, on devrait simplement observer un offset général sur les plots (un offset différent pour chaque ligne, pas de dépendance entre les pas de temps)
SNR=100, SW=1, Δti (zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) + mesure δt avec méthode rapport des phases)
Ref=N (comme ci dessus) Ref=N-1
Dt SNR100 SW1.png
Dt SNR100 SW1 refNm1.png
On a donc simplement appliqué un offset (qui est la différence entre SN et SN-1) à chacune des lignes (à chaque pas de temps)
  • Si on prend successivement toutes les stations comme référence et qu'on fait une moyenne et un écart-type, qu'est ce que ca représente ? Moyenne inférieure aux caractéristiques électroniques des appareils
SNR=100, SW=1, Δti (zoom sur SW avec fenêtre SW synthétiques de largeur variable (uniquement = f(dist,V=[1.5,5]km/s)) + mesure δt avec méthode rapport des phases)
moyenne des Δt en prenant successivement toutes les N stations comme référence écart-type
Dt SNR100 SW1 meanref.png
Dt SNR100 SW1 stdref.png
On a moyenné les Δt obtenus en prenant successivement chacune des stations pour référence. Donc pour chacune de ces station, la dispersion autour de cette moyenne est la même.
Δti si l'on prend successivement les 79 stations comme référence (moyenne et écart-type) (méthode filter pour random noise ; mesure δt avec méthode rapport des phases)
moyenne des Δt en prenant successivement toutes les N stations comme référence écart-type distributions des Δt (avec indication des distributions si l'on ne prend qu'une seule station pour référence = {79,78,1})
SNR=100
Dt SNR100 SW2 meanref.png
Dt SNR100 SW2 stdref.png
Dt gaussianfit SNR 100 SW2.png
SNR=10
Dt SNR10 SW2 meanref.png
Dt SNR10 SW2 stdref.png
Dt gaussianfit SNR 10 SW2.png
SNR=1
Dt SNR1 SW2 meanref.png
Dt SNR1 SW2 stdref.png
Dt gaussianfit SNR 1 SW2.png
Δti si l'on prend successivement les 79 stations comme référence (moyenne et écart-type) (méthode filter pour random noise ; mesure δt avec méthode rapport des phases) + decGF=1
moyenne des Δt en prenant successivement toutes les N stations comme référence écart-type distributions des Δt (avec indication des distributions si l'on ne prend qu'une seule station pour référence = {79,78,1})
SNR=100
Dt SNR100 SW2 decGF1 meanref.png
Dt SNR100 SW2 decGF1 stdref.png
Dt gaussianfit SNR 100 SW2 decGF1.png
SNR=10
Dt SNR10 SW2 decGF1 meanref.png
Dt SNR10 SW2 decGF1 stdref.png
Dt gaussianfit SNR 10 SW2 decGF1.png
SNR=1
Dt SNR1 SW2 decGF1 meanref.png
Dt SNR1 SW2 decGF1 stdref.png
Dt gaussianfit SNR 1 SW2 decGF1.png
  • Si on décale A par rapport à toutes les autres, la solution semble évidente au regard de la matrice antisymétrique des δt, on aura ΔtA décalé et Δti≠A inchangés. De même si on décale A et B. Par contre, si on décale la station choisie comme référence, on biaise tout le monde...
SNR=100 , SW=1 , Experiences de décalage imposé sur les δt(i,j). N est la référence.
δt(1,i)+=1 et δt(i,1)-=1 δt(1,i)+=1 et δt(i,1)-=1 ET δt(2,i)+=-2 et δt(i,2)-=-2 δt(N,i)+=-1 et δt(i,N)-=-1 δt(1,i)+=1 et δt(i,1)-=1 ET δt(2,i)+=-2 et δt(i,2)-=-2 ET δt(N,i)+=-1 et δt(i,N)-=-1
Dt SNR100 SW1 decalage sta1p1.png
Dt SNR100 SW1 decalage sta1p1-sta2m2.png
Dt SNR100 SW1 decalage staNm1.png
Dt SNR100 SW1 decalage sta1p1-sta2m2-staNm1.png
Si la référence est une station saine (N), on retrouve bien Δt(1)=+1 et Δt(i≠1)=0. Si on moyenne les résultats issus de toutes les référence, on introduit une petite erreur (à cause du résultat issu de la référence 1). Ici aussi, on retrouve bien Δt(1)=+1 et Δt(2)=-2 et Δt(i≠{1,2})=0. Ici aussi, en moyennant toutes les références, on introduit une petite erreur. On a imposé un décalage à la station de référence: toutes les autres sont switchées de -décalage. Dans ce cas, si on moyenne les résultats issus de toutes les référence, on a un résultat "juste". On a imposé un décalage à plusieurs stations dont celle de référence. En fait on a un résultat "juste" au décalage de ref près.


SW=2 , Experiences de décalage imposé sur les δt(i,j). N est la référence.
δt(1,i)+=1 et δt(i,1)-=1 δt(1,i)+=1 et δt(i,1)-=1 ET δt(2,i)+=-2 et δt(i,2)-=-2 δt(N,i)+=-1 et δt(i,N)-=-1 δt(1,i)+=1 et δt(i,1)-=1 ET δt(2,i)+=-2 et δt(i,2)-=-2 ET δt(N,i)+=-1 et δt(i,N)-=-1
SNR=100
Dt SNR100 SW2 decalage sta1p1.png
Dt SNR100 SW2 decalage sta1p1-sta2m2.png
Dt SNR100 SW2 decalage staNm1.png
Dt SNR100 SW2 decalage sta1p1-sta2m2-staNm1.png
SNR=10
Dt SNR10 SW2 decalage sta1p1.png
Dt SNR10 SW2 decalage sta1p1-sta2m2.png
Dt SNR10 SW2 decalage staNm1.png
Dt SNR10 SW2 decalage sta1p1-sta2m2-staNm1.png
SNR=1
Dt SNR1 SW2 decalage sta1p1.png
Dt SNR1 SW2 decalage sta1p1-sta2m2.png
Dt SNR1 SW2 decalage staNm1.png
Dt SNR1 SW2 decalage sta1p1-sta2m2-staNm1.png



  • Pour arriver au final à des Δt "absolus", ne pourrait on pas, au final, éliminer des données la station (choisie arbitrairement) de référence. Les N-1 stations restantes seraient considérées comme parfaitement synchronisées ?
  • Je n'arrive toujours pas à comprendre, pourquoi, puisqu'on fait l'hypothèse ΔtN=0, n'aurait on pas le droit de déduire tout bêtement les Δti≠N des δti,N. Dans ce cas on n'utiliserait jamais l'information fournie par les δti,j;i≠N;j≠N.

[edit] 28/03/2018 : Réunion 10h IMAG sur l'étape de prétraitement

Participants : P.-A. Bouttier, A. Lecointre

Durée : 45min

Compte-Rendu :

1) Mémoire à l'exécution du code de prétraitement:

Investigation pb code obspy interpolate qui explose la mémoire:

Voila un bout de code qui génère une trace (remplie de zero) d'une journée à 125Hz et qui interpole linerairement à 100Hz :

  • npts_initial = 86400*125 = 10,800,000 donc 10,800,000*8 = 82MB

Ici on n'a pas besoin de changer le point de départ (minuit pile dans les 2 cas). Le pb de mémoire ne vient donc pas de là.

lecointa@ist-156-56:~/obspy_stream_interpolate$ cat obspy_stream_interpolate.py 
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
from obspy.core import UTCDateTime,Trace
import numpy as np
previous_freq = 125
updated_freq = 100
net = "NE"
sta = "STAT"
loc = "00"
channel = "CHA"
datatype = np.float64
meta = {'station': sta, 'location': loc, 'network': net, 'channel': channel}
mytrace = Trace(data=np.zeros(86400*previous_freq,dtype=datatype),header=meta)
mytrace.stats.sampling_rate = previous_freq
mytrace.stats.starttime = UTCDateTime("2018001")
print(mytrace)
print(mytrace.stats)
nts = mytrace.stats.starttime
#mytrace.interpolate(sampling_rate=updated_freq,starttime=nts,method='linear')
print(mytrace)
print(mytrace.stats)

Quand on commente la ligne interpolate :

SUR LUKE
lecointre@luke:~/IWORMS/IWORMS_preprocess_src$ cat obspy_stream_interpolate/OAR.iworms_memory.5038825.stderr
   81648 Maximum resident set size (KB)
SUR LAPTOP
lecointa@ist-156-56:/tmp/obspy_stream_interpolate$ /usr/bin/time -f "\t%M Maximum resident set size (KB)" python2.7 obspy_stream_interpolate.py >/dev/null
   95024 Maximum resident set size (KB)
SUR ISTOAR
lecoinal@ist-oar:~/tmp/obspy_stream_interpolate$ /usr/bin/time -f "\t%M Maximum resident set size (KB)" python2.7 obspy_stream_interpolate.py > /dev/null
   126324 Maximum resident set size (KB)

Pour interpoler linéairement 125Hz -> 100Hz, 3 working tableaux sont nécessaires. On s'attend donc à consommer 3*80MB = 240MB. Or l'exécution du code consomme 600MB.

Et quand on effectue l'interpolation 125 -> 100Hz :

SUR LUKE
lecointre@luke:~/IWORMS/IWORMS_preprocess_src$ cat obspy_stream_interpolate/OAR.iworms_memory.5038824.stderr
   575164 Maximum resident set size (KB)
SUR LAPTOP
lecointa@ist-156-56:/tmp/obspy_stream_interpolate$ /usr/bin/time -f "\t%M Maximum resident set size (KB)" python2.7 obspy_stream_interpolate.py >/dev/null
   589492 Maximum resident set size (KB)
SUR ISTOAR
lecoinal@ist-oar:~/tmp/obspy_stream_interpolate$ /usr/bin/time -f "\t%M Maximum resident set size (KB)" python2.7 obspy_stream_interpolate.py > /dev/null
   616176 Maximum resident set size (KB)


TODO (P.-A.):

  • Pb mémoire identifié avec toujours obspy-1.0.3 (sur luke env v2 grille avec python/2.7.12_gcc-4.6.2, sur ISTclusters avec python/2.7.13_gcc-6.3.0, sur laptop A.L. avec python/2.7.12_gcc-5.4.0). Comparer avec dernière version stable obspy qui est la 1.1.0 (même si les changelogs ne mentionnent rien en rapport avec notre problème)
  • Creuser dans le code obspy pour identifier la source du problème
  • Faire remonter aux développeurs obspy

2) Parallélisation multicoeurs joblib (et multinoeuds dask ?) du code de prétraitement

Tests de montée en charge WGET SUMMER : Service RESIF@ISTerre m'a relancée et serait très intéressé par la mise en place (rapide) de ces tests (malcommode à effectuer via la grille=contrôle du nb de jobs simultanés). On propose de gérer le parallélisme (sur les traces journalières) au niveau du code python plutôt qu'au niveau du batch.

  • cela impose d'utiliser l'option obspy.fdsn.client pour requêter les données. On laisse donc (temporairement ?) de côté les 3 autres options (wget method GET, wget method POST, accès NFS /media/resif)
    • Note: en séquentiel, on avait mesuré de meilleures perfs avec la wget méthod GET ([[1]]) surtout dans le cas ou il n'y a pas de données. Voir ce que ca donne avec la version parallèle...
  • on commence par juste joblib (mononoeud) sur luke4 ou luke24 ou luke25 (qui sont actuellement les trois noeuds autorisés à effectuer plus de 3 requêtes simultanées)
  • ensuite on ajoutera dask (multinoeud) pour déployer sur luke4+luke24+luke25 (pour augmenter le nb de requêtes simultanées)
  • on va d'abord débrancher les calculs (prétraitement) et l'écriture des outputs HDF5, pour ne tester que les requêtes.
    • Note: quand on voudra rebrancher l'écriture des outputs HDF5, cela va se faire en parallèle : plusieurs "threads" python écriront simultanément, chacun dans un dataset distinct, mais dans le même container HDF5. OK en mononeoud ? OK aussi en multinoeud (je n'ai jamais testé) ?
  • on pourra plus tard se poser la question d'optimiser certaines choses (exemple les coefficients pour les filtres qui pourraient être partagées entre toutes les traces), mais cela risque de nous obliger à quitter obspy ...
  • A priori P.-A. confirme bien que joblib n'imposera pas de synchronisation entre chaque étape de prétraitement. Donc un coeur (=1thread) prend une trace du sac de N traces, la requête, la traite en entier, l'écrit dans un dataset hdf5, et seulement après prend une autre trace du sac de N traces.

TODO:

  • dépôt git (gitlab GRICAD) de la dernière version du code de preprocess (avec le check remove trace si < 2 samples)
  • Fournir un day avec une liste de paramètres (sta,net,chan,loc): en fait c'est fait dans la partie get_metadata du script de soumission actuel

[edit] 14/05/2018 : Réunion 14h ISTerre

Participants : Philippe, Albanne

Durée : 1h30

Ordre du Jour : Préparation de la réunion de vendredi 18/05

ToDo : Pour vendredi, figures à compléter/ajouter dans le CR de la réunion du 26/03 (à déplacer dans une section à part?):

  • Mesure δt: Tester une largeur de fenêtre fixe pour le zoom sur l'onde de surface: on propose de choisir une largeur = 4 x gauspuls_width = 4 x 20sec. En fait j'ai pris 40sec donc 401 points (gauspuls(time,fc=0.15Hz,bw=0.5))...
    • refaires figs d'exemples détaillés ipair=1651 (distance=0km), ipair=1 (distance 30km) et pour ipair=525 (374km)
    • refaire figs Evolution du δt et de la cohérence associée au cours des 24 pas de temps : vérifier que la dépendance en distance disparaît
  • Questionner Laurent Stehly sur ces choix Vmin=1.5km/s, Vmax=5km/sec, window_width=5/f=0.5sec (6 points) ?
  • Distribution des δt mesurés: est ce que l'on fitte mieux la gaussienne si on retire le bin δt=0 ?
  • Une simulation plus réaliste avec décroissance du MaxGF en 1/sqrt(r/r0), au dela d'un certain seuil de distance = 3km/s / 0.15Hz = 20km. Pour r0 on prendrait λ=20km. Ainsi le facteur de décroissance serait:
    • à 20km: 1/sqrt(r/r0)=1/sqrt(20/20)=1
    • à 400km: 1/sqrt(r/r0)=1/sqrt(400/20)=0.22
Decrease synthetic GF amplitude with distance
Decrease GF with distance.png
  • Δt:
    • dupliquer toutes les figures pour SNR=10 et SNR=1 (y compris les figs avec décalae imposé)
    • ajouter figs distribution des Δt dans le cas "une station de référence" (en choisir qqunes) et dans le cas ou l'on équilibre toutes les stations (cas "toutes les stations choisies successivement comme référence"): vérifier que la distribution sera mieux randomisée et plus étroite.
  • Complétude: essayer de déterminer à la main le ou les sous-graphe sur un cas idéal avec faible nombre de stations et ou l'on efface volontairement une ou des arêtes, déterminer à la main (sur la matrice antisymétrique) le ou les sous-graphes restants...


Discussions:

  • Mesure de différences d'horloge inter-stations: δt+ et δt- : séparer parties physique ((δt+ + δt-)/2) et non physique ((δt+ - δt-)/2), et faire les 2 calculs Δt+ et Δt- associés.
  • Retour aux données réelles: est ce que l'on abandonnerait (temporairement) les données AlpArray pour revenir aux données Pyrope (déjà étudiées par Laurent). Dans ce cas, référence periode = 1 an et timestep = 1 jour.

[edit] 18/05/2018 : Réunion 13h 158 IMAG

Participants : G. Arneodo, J. Schaeffer, C. Picard, V. Louvet, P.-A. Bouttier, P. Roux, A. Lecointre

Durée : 2h30

Compte-rendu :

1) 30' : Prétraitements : déploiement joblib (?) et monitoring accès Resif@Summer.

Invités ISTerre-GeoData-RESIF: G. Arneodo (resp. Exploitation et production, dev. logiciels monitoring, SGBD), J. Schaeffer (resp. Infrastructures matérielles et middleware centre de données RESIF)

[support de discussion du 28/03 (P.-A.B., A.L.)]

Deux points à tester:

  • monitoring de l'"accès réel aux données": être capable de faire un retour (quelles données ai je pu télécharger ?) versus (quelles données étaient disponibles?). Actuellement les logs d'exécution du code de pretraitement indiquent, jour par jour, et trace par trace, si au moins un échantillon (pas forcément la journée entière) a bien été téléchargé. On ne sait pas si la donnée est partielle ou complète. Exemple:
lecointre@luke:~/IWORMS/IWORMS_preprocess_grid$ cat OAR.cigri.10970.3759238.stdout
2017056
//bettik/lecointre/IWORMS_preprocess_10970_3759238_2017056
RC du wget: 200
wget : 1 query for each daytrace
FR ARBF 00 HHZ : 1 / 73 : 15.223298375 sec
FR ARTF 00 HHZ : 2 / 73 : 12.065697019 sec
FR BANN 00 HHZ : 3 / 73 : 8.087988037 sec
FR BLAF 00 HHZ : 4 / 73 : 9.906879376 sec
...
FR SURF 00 HHZ : 44 / 73 : 10.012187251 sec
FR TRBF 00 HHZ : 45 / 73 : Unknown format for file /tmp/obspy-8IzVuK.tmp
4.652756058 sec
FR TRIGF 00 HHZ : 46 / 73 : 7.971614358 sec
...

Sur cet exemple le log indique que FR.TRBF.00.HHZ n'a pas délivré de données entre 2017055.23:59:00.00 et 2017057.00:01:00.00, cette station sera absente de notre mesure de δt pour cette journée 2017056. Le log indique également que les autres stations ont délivré de la donnée partielle (plusieurs morceaux séparés par des trous: dans ce cas le code de prétraitement bouche les trous) ou complète, ces stations sont donc utilisées pour la mesure du δt.

Je suppose qu'il faudrait donc que les logs indiquent, trace par trace et jour par jour, entre quelles dates précises les données ont été téléchargées (= quels trous ont été bouchés) ? Sous quel format l'indiquer pour être exploitable facilement par les gens RESIF pour contrôle ?

  • monitoring de la charge RESIF@SUMMER occasionné par le déploiement grande échelle de ce prétraitement: pour cela on a besoin d'une solution de parallélisation autre que la grille CIMENT (afin de pouvoir contrôler parfaitement le nombre de processus simultanés).

Décision: mise en oeuvre de deux solutions comparatives de parallélisation du prétraitement Python:

  • joblib (mononoeud) associé à dask (multinoeud)
  • MPI4py (mononoeud et multinoeud)

P.-A. B. se charge d'avancer sur ces développements d'ici la prochaine réunion prévue en juillet. Pour aller au plus efficace/simple sur ces prétraitements, ne retenir que les parties suivantes du code:

  • récupération des données d'input via l'option obspy.fdsn.client, laisser de côté les 3 autres options (1) wget GET method avec flux "stream", (2) wget GET method avec avec écriture fichier temporaire, (3) wget POST method (3 est non fonctionnelle).
  • directement filtre 1-Bit sur données inputs (on "oublie" les prétraitements précédents)
  • écriture // dans dataset HDF5

Questions:

  • quel support OAR pour dask ? aucun par défaut.
  • les IO // avec joblib/dask ?

-> intérêt de l'approche comparative joblib/dask versus MPI4py.

2) 15' : Prétraitements : mémoire

[support de discussion du 28/03 (P.-A.B., A.L.)]

Après analyse il s'avère que le problème d'empreinte mémoire vient de l'interpolation SciPy v1.1. (mémoire * 8).

Pour info: interp linéaire numpy.interp (mémoire * 6.5).

Une solution serait de switcher vers un filtre Lanczos (recodé dans obspy). Ainsi on ne dépasse pas les 240MB (mémoire * 3 : OK).

[Lanczos resampling (wikipedia)]

[Doc obspy lanczos resampling]

Questions:

  • check eventual resulting NaN or Inf values ? A.L. indique avoir eu précédemment des ennuis avec la méthode weighted_averaged_slopes (SAC), c'était la raison du passage à la méthode linear.
  • quel choix de paramètre pour la taille de la fenêtre pour le sinc ? Une valeur a=20 -> T_CPU * 10.
  • comment ca va se passer quand le segment a très peu de points (<20) ? Actuellement on garde tous les morceaux dès lors que nb_points>=2.

P.R. et A.L. feront des tests pour déterminer a (qui sera peut être donnée-dépendant), et pour mesurer l'impact de ce changement de filtre sur les mesures δt des données réelles/synthétiques.

3) 1h45' : δt et Δt synthétiques

[Support de discussion du 26/03 (C.P., P.R., A.L.)] et [support de discussion du 14/05 (P.R., A.L.)]

Pour conclure cette étude sur les données synthétiques, on veut encore faire des plots statistiques de l'evolution de la moyenne et variabilité quand on modifie:

  • le nombre de stations (actuellement 79)
  • le nombre de pas de temps dans la période de référence (actuellement 24)

+ ajouter un paramètre de décroissance de l'amplitude de la GF avec la distance.

+ décalage appliqué dans les signaux synthétiques plutôt que dans les δt mesurés.

+ les distributions de Δt "paraissent" ordonnées par ordre de station de référence. ca devrait être aléatoire. A vérifier.

Question : on n'avait pas fait de 1bit dans les signaux synthétiques avant mesure δt. A tester aussi ? ca va impacter le contrôle du SNR ?

Rq: il doit y avoir une coquille "distance" versus "distance/2" -> pas d'impact sur notre analyse mais à corriger.

[edit] 07/06/2018 : Réunion-Discussion slides de Chamonix dans le cas graphe non complet

Participants : C. Picard, A. Lecointre

Durée : max 1h

Compte-rendu :

1) Cas graphe complet. Globalement tout est OK. Convention Nb noeuds = n+1.

  • slide 5 : GT G m = GT d
  • slide 6 : v0 = 1/(n+1) J1,n+1 pour avoir norme(v0)=1
  • slide 7 : on passe bien dans le RHS VT et pas V-1 (intérêt de ne pas avoir à inverser). A noter que VT n'est cependant pas rigoureusement égal à V-1 (première ligne avec des 1/(n+1) vs première ligne avec des 1). Ca n'a pas d'importance car ca disparaît avec la première valeur propre λ0=0.
  • slide 8 : m(1)=α (α=0) et D(k) (VT m)(k) = b(k) = (VT GT d)(k)
  • slide 9 : les VP sont 0 et 4 (si n+1=4noeuds)

2) Cas graphe non complet.

  • slide 12 : Rq : Quelquesoit la complétude du graphe, la matrice des degrés s'obtient en faisant la somme des éléments de chaque ligne (ou colonne si non orienté) de la matrice d'adjacence. Or la matrice d'adjacence se construit très facilement à partir de notre base d'observations : c'est bool(δti,j). Ca va être beaucoup plus intéressant de construire la Laplacien ainsi plutôt qu'à partir de la matrice d'incidence du graphe G.
  • slide 13 : Définition: graphe connexe = on n'a pas plusieurs graphes disjoints.
  • slide 14. I aurait pour dimension [4,12=2*6] si 4 noeuds (car 12 arêtes orientées possibles avec 4 noeuds). I contient des 0 (=pas d'arête) et des -1 (=vers). I contient 4 blocs : sort de 1, sort de 2, sort de 3, sort de 4. Dans le cas d'un graphe complet orienté avec 4 noeud, chaque noeud est 3 fois un puit. Pb car (si gr complet) IT I ne redonne pas (D-A) le laplacien du graphe complet : on ne retrouve pas cette autre définition du laplacien pour le cas d'un graphe orienté : on retrouve bien les 3 sur la diagonale mais on a perdu tous les -1 (deviennent 0)...
  • slide 16 : Rq : cas d'un graphe orienté non complet dans lequel il manque l'arête de 1 vers 2. La matrice des degrés devient diag(3,2,3,3). cette matrice des degrés seule ne renseigne pas sur l'arête manquante, on sait juste que c'est une des 3 arêtes qui devraient arriver en 2. C'est la matrice d'adjacence qui va fournir cette info. De plus dans le cas du graphe orienté la matrice d'adjacence n'est plus nécessairement symétrique (on peut avoir existence d'une arête dans un sens mais pas dans l'autre).
  • on ne va pas au delà de la slide 16 tant qu'on n'a pas clarifié l'expression de I. Fin de la réunion.

[edit] 28/06/2018 : Réunion Logging prétraitements 10h 245 ISTerre

Participants : G. Cougoulat, G. Arneodo, A. Lecointre

Durée : 1h

Compte-Rendu :

  • les bonnes pratiques de logging des prétraitements iworms pour l'exploitation des logs pour le monitoring de l'accès aux données RESIF@SUMMER (mots clés : Elasticsearch , Zabbix , Kibana)

Présentation de Gregory: File:Logging@resif.pdf

Actuellement:

  • logs "exportés" sur serveurs distants (sécurité)
  • plateforme de test sur Winter
    • Elasticsearch : DB
    • Logstash: logiciel pour alimenter DB à partir de logs
    • Kibana (Graphana) : logiciel de visu
  • une VM Winter avec 250GB espace disk

Autres besoins émergent : DashBoard (responsable scientifiques et techniques RESIF), Colmet (monitoring de jobs OAR), ...

Graylog2 (fork ELK) : opensource, mais non utilisés par OSUG ni par collègues européens (EIDA)

Contrainte de compatibilité avec EIDA (collègues Europe RESIF).

WFCatalog, code Python qui déroule toute la chaîne, en cours de test.

librairie logger de python

Dans le cadre de nos futurs déploiements luke de rapatrie+preprocess de données RESIF@SUMMER (dans le cadre du projet iWORMS mais pas que), pour exploiter ces campagnes de stress-test, il serait souhaitable de sortir dans les "logs" (dans les stdout dans un premier temps):

  • date de début de job OAR
  • date de fin de job OAR
  • url de la requête wget
  • nb d'échantillons récupérés
  • taille du fichier MSEED récupéré: à noter qu'à priori on ne dispose pas de cette info puisque notre workflow est :
# option wget
wget $url -O - | python read_and_preprocess.py
# ou option obspy.fdsn.client
python rapatrie_and_read_and_preprocess.py

Dans les deux cas la donnée requêtée n'est jamais écrite sur disque, le flux est directement chargé en mémoire et prétraité. On prévoit donc de pouvoir jouer ce test aussi (même si ca ne sera pas utilisé en production):

wget $url -O - > tmp.mseed
ls -l tmp.mseed
cat tmp.mseed | python read_and_preprocess.py

A noter que je ne sais pas trop comment mener nos campagnes de stress-test sans pouvoir vider les cache RESIF ET les cache luke: exemple aujourd'hui j'observe:

wget $url -O - > /dev/null  # 5sec elapsed time
# on relance ultérieurement (même si jobs OAR différents)
wget $url -O - > /dev/null  # 0.3sec elapsed time

Si on veut jouer la campagne de stress-test sur qq 100aines de coeurs luke, il faudra prévoir autant de requêtes url différentes, et prévoir un temps d'attente avant d'éventuellement rejouer le même stress-test.

A noter que Gregory peut vider les cache RESIF. Par contre on ne peut pas vider les cache luke à la demande (?)

Actions concrètes à mettre en oeuvre à l'issue de la réunion:

  • Albanne fournit à Gregory un exemple de stdout avec la liste de champs ci dessus, écrits "en vrac" à la "echo", juste pour un seul accès (une seule url requête)
  • Gregory reprend le fichier "à la main" et le réécrit dans un format compatible syslog "formatter_EidaFormatter" (pour fournir un modèle).
  • Albanne, soit reprend les instruction "echo" du script de soumission, soit fabrique un script shell pour traduire automatiquement (afin de ne pas casser les outils actuels à la gnuplot pour exploiter les infos de perfs des stdout).
  • Voir avec P.-A. pour l'avancement de la parallélisation MPI4Py + joblib/dask du code de preprocess (du moins juste les entrées sorties de ce code pour faire au plus rapide/facile), afin d'avoir un contrôle strict sur le nombre d'instructions python rapatrie_and_read_and_preprocess.py (contrôle qu'on n'a pas actuellement avec nos déploiements CiGri)

Rq : la semaine dernière dans le cadre des travaux autour du filtre Lanczos, j'ai redéployé le prétraitement pour un an de données (~70 stations 43N-46N-3E-8E) sur la grille (mais sans le max_jobs=3 comme je l'avais fait l'an dernier, puisque maintenant la limitation à 3 requêtes simultanées a été levée), et j'avais quelques 365 jobs cigri simultanés (un job cigri par paramètre:($YEAR,$JDAY)) qui faisaient:

for urltrace in $listurltrace; do
  python rapatrie_and_read_and_preprocess.py $urltrace $YEAR $JDAY
done

et ca a tenu.

Exemple de stdout pour Gregory :

1) cas interp linear, rapatrie option = obspy.fdsn.Client (pas de fichier tmp.mseed sur disque) :

$ cat OAR.logging.7602869.stdout
luke4
2017011
//bettik/lecointre/IWORMS_preprocess_00000_7602869_2017011
RC du wget: 200
obspy.clients.fdsn.Client
FR ARBF 00 HHE : 1 / 3 :
start 1530187935.589382686
('client.get_waveforms', 'FR', 'ARBF', '00', 'HHE', UTCDateTime(2017, 1, 10, 23, 59), UTCDateTime(2017, 1, 12, 0, 1))
('total_npts', 10815412)
end 1530187949.912103892
elapsed time 14.322721206 sec
2017 011 FR ARBF 00 HHE
FR ARBF 00 HHN : 2 / 3 :
start 1530187950.085340237
('client.get_waveforms', 'FR', 'ARBF', '00', 'HHN', UTCDateTime(2017, 1, 10, 23, 59), UTCDateTime(2017, 1, 12, 0, 1))
('total_npts', 10815852)
end 1530187965.636257654
elapsed time 15.550917417 sec
2017 011 FR ARBF 00 HHN
FR ARBF 00 HHZ : 3 / 3 :
start 1530187965.726702849
('client.get_waveforms', 'FR', 'ARBF', '00', 'HHZ', UTCDateTime(2017, 1, 10, 23, 59), UTCDateTime(2017, 1, 12, 0, 1))
('total_npts', 10816366)
end 1530187981.206002101
elapsed time 15.479299252 sec
2017 011 FR ARBF 00 HHZ

2) cas interp linear, rapatrie option = wget $url -O - > tmp.mseed ; ls -l tmp.mseed ; cat tmp.mseed | python read_and_preprocess.py :

$ cat OAR.logging.7602871.stdout
luke4
2017011
//bettik/lecointre/IWORMS_preprocess_00000_7602871_2017011
RC du wget: 200
wget : 1 query for each daytrace
start 1530188018.899082162 
FR ARBF 00 HHE : 1 / 3 : 
-rw-r--r-- 1 lecointre l-isterre 23474176 Jun 28  2018 tmp.mseed
('total_npts', 10815412)
end 1530188028.126085170 
elapsed time 9.227003008 sec
url: http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHE&starttime=2017-01-10T23:59:00.000000&endtime=2017-01-12T00:01:00.000000
start 1530188028.220952213
FR ARBF 00 HHN : 2 / 3 :
-rw-r--r-- 1 lecointre l-isterre 23494656 Jun 28  2018 tmp.mseed
('total_npts', 10815852)
end 1530188037.320097179
elapsed time 9.099144966 sec
url: http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHN&starttime=2017-01-10T23:59:00.000000&endtime=2017-01-12T00:01:00.000000
start 1530188037.407966422
FR ARBF 00 HHZ : 3 / 3 :
-rw-r--r-- 1 lecointre l-isterre 23355392 Jun 28  2018 tmp.mseed
('total_npts', 10816366)
end 1530188046.020728761
elapsed time 8.612762339 sec
url: http://ws.resif.fr/fdsnws/dataselect/1/queryauth?network=FR&station=ARBF&location=00&channel=HHZ&starttime=2017-01-10T23:59:00.000000&endtime=2017-01-12T00:01:00.000000

Les "echo date" sont au format $(date +%s.%N) : pratique pour faire du bc , mais possible de changer.

[edit] 06/07/2018 : Réunion 9h30 ISTerre pour préparer réunion du 10/7

Participants : P.Roux, A.Lecointre

Durée : 1h

Compte Rendu : Préparation réunion du 10/07 (voir support de discussion ci dessous)

  • Améliorer les plots : fit gaussien -> trouver une représentation de uniquement (on sait que moyenne est nulle) σ en fonction de SNR, nstep, Nb_stations, (decGF, ) ...
  • Cas nbpair=27966 (237 "stations": en fait 79 mais avec les 3 composantes):
    • vérifier quelle largeur SW quand distance=0 (40sec ou 20sec ?)
    • nbpair = 3081 ou 27966: actuellement 2 familles de problèmes différents (augmentation de densité). Peut être plutôt reprendre une box de grille de données réelles augmentée de sqrt(3) au nord et à l'est (centre France), pour espérer avoir 3 fois plus de stations toutes HHZ sans augmenter artificiellement la densité. Ou alors (pas dans l'immédiat) créer une grille synthétique adaptable (distance, densité) à souhait.
  • Pour début septembre, préparer quelques plots pour habiller un article sur: données synthétiques, statistiques, graphe complet, erreurs de temps synthétiques. garder en ouverture le cas graphe incomplet et les données réelles (SanJacinto?)
  • Pour début septembre: revenir à un cas pratique, on abandonne AlpArray, on utilise San jacinto (1108 Z nodes covering <1km2, instrument spacing of 10-30 m). Box ~100 stations au coeur du réseau. Viser jday=150 et 2-4Hz : c'est à dire expérience A03 ([2]) où on a déjà les GF 10min et moyenne journalière (donc nstep=144); Prévoir éventuellement d'étendre A03 sur quelques julianday avant et après (demander à Chloé sur quelles dates a t on le max de stations en fonctionnement ?, ou le déduire de mes logs de prétraitement : 1) [pour identifier les JDAY en ilieu d'expérience avec beaucoup de gaps] et 2 [pour éliminer les JDAY avec peu de stations: par ex ne conserver que 130<=jday<=154]
  • Expériences de décalage imposé: comprendre pourquoi ca ne semble pas fonctionner à l'augmentation du nb de station ou de pas de temps (les toutes dernières figures). pour le test 2 et 400 pas de temps, faire le décalage variable 0.1sec par pas de temps uniquement pendant les 24 premier pas de temps. ne pas représenter les pas de temps ou il n'y a pas de décalage (moins surcharger les plots). test avec SNR=100 pour validation.

[edit] 10/07/2018 : Réunion 14h en 258 IMAG

Participants : P.-A. Bouttier, C.Picard, A.Lecointre

Durée : 45min

Compte-Rendu :

  • 1.1)
    • Fournir des données 1bit ainsi que des dt et coherence pour des expériences données réelles et synthétiques pour l'écode d'été de Christophe. Données et infos ici : luke:/bettik/lecointre/IWORMS/WORKSHOP_SERBIE
  • 1.2)
    • Mise à jour du dépôt Git du code de prétraitement (avec le print des infos logging pour Gregory)
    • Clarifier avec Gregory ce besoin(?) de l'info de taille en MB du fichier MSEED téléchargé : incompatible avec notre workflow.
    • P.-A. prépare une version MPI4Py pour la prochaine réunion mi-fin septembre
  • 2)
    • Constat intéressant de la dégradation des perfs d'accès RESIF@SUMMER
    • Faire ticket à RESIF pour signaler ce problème de flag quality='M' dont les effets diffèrent entre wget et obspy
    • Pousser les tests Lanczos a=5 voire moins (avec contrôle nbdif_1bit / interp linéaire)
  • 3)
    • Préparer "plots synthèse" de σSNR*Δt = f(SNR,nbStations,densiteStation,nstep,decGF) pour illustrer papier méthodologie
    • Essayer de mesurer et représenter l'incertitude sur Δt en fonction de (SNR,nbStations,densiteStation,nstep,decGF)
    • Mise en évidence de cas ou la méthode peut donner un résultat faux malgré un très bon niveau de cohérence : quand l'une des station est désynchronisée pendant toute ou une trop grande partie de la période de référence (voir figures exemple dans l'OdJ).

Ordre du jour :

1) Point d'info sur les travaux qui ont eu lieu en comité restreint depuis dernière réunion:

  • 1.1) Graphes complets, non complets [07/06/2018]
  • 1.2) Logging des téléchargements et prétraitements [28/06/2018]

2) Prétraitements : Interpolation -> 100Hz : pb conso mémoire interp linéaire. Test interp Lanczos.

Test 2017011, données réelles 75 stations dont 71 avec des données. Box 43N-46N-3E-8E.

RAPATRIE_OPTION=1
JJJ="2017011"
listnet="FR,RD,RA,ZH,Z3,G"
liststa="*"
listloc="*"
listchan="HHZ"
minlat=43
maxlat=46
minlon=3
maxlon=8

-> ne devrait on pas plutôt tester Lanczos sur des données synthétiques ?

On fixe les paramètres ci dessus, on fixe la méthode RAPATRIE_OPTION = 1 (obspy.fdsn.Client)

2.0) Le pb initial de consommation mémoire semble résolu : passage de 650MB à 400MB (pour de la données ~100Hz journalière)

2.1) Imact sur les données : on mesure le pourcentage de variations sur les données 1-Bit :

  • Lanczos/Linear : un max 0.018%, en général de l'ordre de 1e-4%

2.1) on rejoue la méthode linear : dégradation perfs 3-6sec (plot de droite croix bleues) -> 10sec (plot de gauche croix bleues)

[Ticket resif-dc #0106441]

Réponse proposée : la limitation du nombre de connection en vigueur en 2017 a été levée, la machine ws.resif.fr est beaucoup plus chargée maintenant. Il y a bien plus d'accès à SUMMER. Il y a également plus de connexions parallèles au SGBD.

2.2) lanczos avec a = 20 : 30 ou 150 sec / trace

2.3) lanczos avec a = 10 : 20 ou 80 sec / trace

Pourquoi certaines traces coutent plus cher que d'autres à prétraiter avec méthode lanczos ? On n'avait pas fait cette constatation avec méthode linear. Ca n'est pas lié à la présence de trous ou au nombre de points à traiter.

2.4) téléchargement seul (wget) : 5-6sec (wget -O - > /dev/null : points bleu cyan) si la donnée n'est ni en cache luke ni en cache RESIF. -> A retenir qu'il faut mettre une barre d'erreur de 5-6sec autour de tous ces points.

2.5) téléchargement + read (qqs wget ou obspy.fdsn.Client) : 11-12 sec (ronds noirs et points verts) : avec donnée en cache luke et en cache RESIF

2.6) On joue le rapatrie+preprocess lanczos sur un an de données:

rappel: Get Input + Preprocess 2017011, 216 traces
Lanczos linear elps.png
PreprocessIWORMS NFS WS.png
Lanczos linear resmem.png
Lanczos linear npts.png
Quality.png
Lanczos linear dif.png

3) δt et Δt synthétiques

3.1) on complète les plots discutés lors de la réunion du 18/05/2018 (plots dans le CR de la réunion du 26/03/2018) : Δti si l'on prend successivement les 79 stations comme référence (moyenne et écart-type) (méthode filter pour random noise ; mesure δt avec méthode rapport des phases) et decGF=0/1 :

  • on veut vérifier cette suspicion de gaussienne ordonnée selon numéro de station
  • on veut voir l'effet de la decroissance GF en fonction de la distance (decGF=0/1):

Rappel: decGF=1 : Une simulation plus réaliste avec décroissance du MaxGF en 1/sqrt(r/r0), au dela d'un certain seuil de distance = 3km/s / 0.15Hz = 20km. Pour r0 on prendrait λ=20km. Ainsi le facteur de décroissance serait:

  • à 20km: 1/sqrt(r/r0)=1/sqrt(20/20)=1
  • à 400km: 1/sqrt(r/r0)=1/sqrt(400/20)=0.22
Decrease synthetic GF amplitude with distance
Decrease GF with distance.png


Δti si l'on prend successivement les 79 stations comme référence [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2)]
decGF=0 decGF=0 decGF=1 decGF=1
SNR=100
Dt allref gaussianfit SNR 100 SW2 decGF0.png
Dt allref gaussianfit3D SNR 100 SW2 decGF0.png
Dt allref gaussianfit SNR 100 SW2 decGF1.png
Dt allref gaussianfit3D SNR 100 SW2 decGF1.png
SNR=10
Dt allref gaussianfit SNR 10 SW2 decGF0.png
Dt allref gaussianfit3D SNR 10 SW2 decGF0.png
Dt allref gaussianfit SNR 10 SW2 decGF1.png
Dt allref gaussianfit3D SNR 10 SW2 decGF1.png
SNR=1
Dt allref gaussianfit SNR 1 SW2 decGF0.png
Dt allref gaussianfit3D SNR 1 SW2 decGF0.png
Dt allref gaussianfit SNR 1 SW2 decGF1.png
Dt allref gaussianfit3D SNR 1 SW2 decGF1.png

Conclusions du 3.1):

  • la forme de la distribution des Δt n'est pas ordonnée par ordre de numéro de station. ouf!
  • quand on fait décroître l'amplitude de la GF synthétique avec la distance, la distribution des Δt s'affaisse (plus large et moins haute) -> on abaisse le rapport signal à bruit.
  • Calculer le Δt en moyennant les Δt obtenus en prenant successivement pour référence chacune des stations permet d'affiner et réhausser la distribution gaussienne des Δt. Cela reste vrai que l'on ait fait décroître ou non l'amplitude de la GF synthétique avec la distance.

3.2) Etude statistique avec les données synthétiques.

On fixe les paramètres:

  • CC=0 : méthode rapport des phases pour mesure δt
  • SW=2 : zoom sur SW avec fenêtre SW synthétiques de largeur fixe 40sec
  • decGF=1 : on applique la décroissance de l'amplitude de la GF synthétique en fonction de la distance.

On joue sur ces trois paramètres:

  • SNR : {1;10;100}
  • nb de pas de temps dans la moyenne de référence: nstep : actuellement 24 ; test 720 (=30*24h=1mois) ; test 1440 (2mois).
  • nbpair : actuellement 3081=79*(79-1)/2 (79 stations only HHZ). test avec les 3 composantes HH[ENZ] = (3*79)*(3*79-1)/2=27966 (on bidouille le fichier distance ainsi: pour les 3081 paires HHZ on garde les mêmes distances mais on duplique 9 fois chaque paire, puis pour les autopaires on rajoute trois fois chaque autopaire en affectant distance=0).

3.2.1) Distribution des δt

nbpair=3081: δti,j [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1]
nstep=24 nstep=720 nstep=1440
SNR=100
Dt distribution SNR 100 SW2 CC0 FILTER decGF1 24 3081.png
Dt distribution SNR 100 SW2 CC0 FILTER decGF1 720 3081.png
Dt distribution SNR 100 SW2 CC0 FILTER decGF1 1440 3081.png
SNR=10
Dt distribution SNR 10 SW2 CC0 FILTER decGF1 24 3081.png
Dt distribution SNR 10 SW2 CC0 FILTER decGF1 720 3081.png
Dt distribution SNR 10 SW2 CC0 FILTER decGF1 1440 3081.png
SNR=1
Dt distribution SNR 1 SW2 CC0 FILTER decGF1 24 3081.png
Dt distribution SNR 1 SW2 CC0 FILTER decGF1 720 3081.png
Dt distribution SNR 1 SW2 CC0 FILTER decGF1 1440 3081.png
nbpair=27966: δti,j [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1]
nstep=24 nstep=720 nstep=1440
SNR=100
Dt distribution SNR 100 SW2 CC0 FILTER decGF1 24 27966.png
Dt distribution SNR 100 SW2 CC0 FILTER decGF1 720 27966.png
Dt distribution SNR 100 SW2 CC0 FILTER decGF1 1440 27966.png
SNR=10
Dt distribution SNR 10 SW2 CC0 FILTER decGF1 24 27966.png
Dt distribution SNR 10 SW2 CC0 FILTER decGF1 720 27966.png
Dt distribution SNR 10 SW2 CC0 FILTER decGF1 1440 27966.png
SNR=1
Dt distribution SNR 1 SW2 CC0 FILTER decGF1 24 27966.png
Dt distribution SNR 1 SW2 CC0 FILTER decGF1 720 27966.png
Dt distribution SNR 1 SW2 CC0 FILTER decGF1 1440 27966.png

Conclusions du 3.2.1):

La distribution des δt suit une loi gaussienne:

  • impact SNR : σ (quand SNR diminue, les niveaux de coherence baissent et l'amplitude du max(δt) augmente)
  • impact nstep : fit des données à la gaussienne (visible sur les faibles SNR)
  • le nb de stations n'a pas d'impact sur les δt (ce qui est tout à fait normal, le calcul du petit δt pour une paire donnée est totalement indépendant du calcul pour une autre paire).

3.2.2) Distribution des Δt

nbpair=3081: Fit gaussien de la distribution des Δti en prenant successivement chacune des N stations pour référence, puis en moyennant ces N résultats [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1]
nstep=24 nstep=720 nstep=1440
SNR=100
Dt gaussianfit3D SNR 100 SW2 decGF1 nstep24 nbpair3081.png
Dt gaussianfit3D SNR 100 SW2 decGF1 nstep720 nbpair3081.png
Dt gaussianfit3D SNR 100 SW2 decGF1 nstep1440 nbpair3081.png
SNR=10
Dt gaussianfit3D SNR 10 SW2 decGF1 nstep24 nbpair3081.png
Dt gaussianfit3D SNR 10 SW2 decGF1 nstep720 nbpair3081.png
Dt gaussianfit3D SNR 10 SW2 decGF1 nstep1440 nbpair3081.png
SNR=1
Dt gaussianfit3D SNR 1 SW2 decGF1 nstep24 nbpair3081.png
Dt gaussianfit3D SNR 1 SW2 decGF1 nstep720 nbpair3081.png
Dt gaussianfit3D SNR 1 SW2 decGF1 nstep1440 nbpair3081.png
nbpair=27966: Fit gaussien de la distribution des Δti en prenant successivement chacune des N stations pour référence, puis en moyennant ces N résultats [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1]
nstep=24 nstep=720 nstep=1440
SNR=100
Dt gaussianfit3D SNR 100 SW2 decGF1 nstep24 nbpair27966.png
Dt gaussianfit3D SNR 100 SW2 decGF1 nstep720 nbpair27966.png
Dt gaussianfit3D SNR 100 SW2 decGF1 nstep1440 nbpair27966.png
SNR=10
Dt gaussianfit3D SNR 10 SW2 decGF1 nstep24 nbpair27966.png
Dt gaussianfit3D SNR 10 SW2 decGF1 nstep720 nbpair27966.png
Dt gaussianfit3D SNR 10 SW2 decGF1 nstep1440 nbpair27966.png
SNR=1
Dt gaussianfit3D SNR 1 SW2 decGF1 nstep24 nbpair27966.png
Dt gaussianfit3D SNR 1 SW2 decGF1 nstep720 nbpair27966.png
Dt gaussianfit3D SNR 1 SW2 decGF1 nstep1440 nbpair27966.png
Distribution et fit gaussien des Δti obtenus en moyennant les N Δti1,i2 provenant du calcul successif en prenant chaque station pour référence [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1]
Dt gaussianfit meanref SW2 decGF1.png
Dt R8 gaussianfit meanref SW2 decGF1.png

Conclusions du 3.2.2):

La distribution des Δt suit une loi gaussienne:

Sur les Δt en prenant une station pour référence:

  • qqsoient les paramètres nstep, SNR, nb_stations, on confirme que:
    • pas d'ordonnancement de la distribution selon le numéro de la station
    • σ diminue quand on moyenne les N références

Sur les Δt obtenus en moyennant les N références:

  • impact du SNR: σ (et donc amplitude de abs(Δt)) est une fonction linéaire du SNR.
  • impact du nstep: la distribution fitte mieux la loi gaussienne
  • impact du nombre de stations à densité ~ constante : non visible sur SNR=100 et 10. σ diminue légèrement sur SNR=1.
  • impact de la densité de stations : σ diminue

Il faudrait envisager une représentation de σSNR*Δt (plutot que σΔt) en fonction de nstep, nbsta ou densité, decGF, ...

3.3) On applique un décalage sur une station dans les données synthétiques (GF) plutôt que directement dans les δt mesurés à partir de données synthétiques.

On fixe:

  • nbpair = 3081 -> 79 stations
  • SW=2 (largeur fixe 40sec = 401 samples)
  • CC=0
  • decGF=1
  • nstep=24

On fait 2 tests:

  • test 1 : on décale la station 1 d'1 sec (10 samples) au 1er pas de temps uniquement
  • test 2 : pendant les 24 premiers pas de temps, on décale la station 1 progressivement à chaque pas de temps de 0.1sec*istep, donc de 0.1sec (1 sample) à 2.4sec (24 samples), puis on maintient ce décalage fixe de 2.4sec jusqu'à la fin de la période de référence
Moyenne des Δt obtenus en prenant successivement chacune des N=79 stations pour référence. nbpair=3081, SW=2, decGF=1, méthode filter pour random noise, mesure δt avec méthode rapport des phases (CC=0)
test 1 test 2
Dt decalage FIXE R8 SW2 decGF1.png
Dt decalage ISTEP R8 SW2 decGF1.png
L'incertitude sur Δt devient inacceptable quand SNR=1. Si une station est désynchronisée pendant toute la période de référence, la méthode donne un résultat faux (pour cette station).


nstep=24 -> 720 et nbsta=79 -> 130:
test 1 test 2
SW=100
Dt decalage FIXE R8 SNR100 SW2 decGF1.png
Dt decalage ISTEP R8 SNR100 SW2 decGF1.png
SW=1
Dt decalage FIXE R8 SNR1 SW2 decGF1.png
Dt decalage ISTEP R8 SNR1 SW2 decGF1.png
L'incertitude sur Δt reste inacceptable quand SNR=1, même si on aumente le nombre de step ou le nombre de stations (à densité ~ constante). Même si on augmente le nombre de stations (à densité ~ constante), la méthode donne une résultat faux (pour cette station) si une station est désynchronisée pendant toute la période de référence.


Conclusions du 3.3):

  • incertitude sur Δt (voir ccl dans tableau ci dessus)
  • cas d'une désynronisation longue période (voir ccl dans tableau ci dessus)

[edit] 01-10-2018 : Reunion de Rentree 9h IMAG 258

Participants : P.-A. Bouttier, C. Picard, P. Roux, L. Stehly, D. Wolinieck, C. Pequegnat, J. Touvier, J. Schaeffer, G. Arneodo, A. Lecointre

Durée : 3h

Compte-Rendu (le point 5 intègre aussi le CR d'un petit debriefing du 02/10/2018 entre P.R., L.S., A.L.) :

1) Point d'info prétraitement interp Lanczos

  • interp linear scipy v1.1 : mem*8 (650MB pour l'ensemble du preprocess) / interp Lanczos a={20,10,5} : mem*3 (400MB)
  • traité 1 an de 80 stations en Lanczos : pas vu de NaN ni Inf
  • impact Lanczos a={20,10,5} sur %dif dans data_1Bit : max 30 à 34 pts diffèrent par rapport à linear sur une trace de 8,640,000 points

Conclusion : on choisit Lanczos avec a=5

2) Logging prétraitements

2 types de requêtes dans le cahier des charges pour la version MPI4Py : "à la wget" et "obspy"

pour le "à la wget" (qui permettrait de jouer avec le champ refresh_cache dans l'url, chose qu'on ne peut pas faire en obspy), proposition d'utiliser urllib, comme cet exemple dans la doc du portal resif :

# this will make a wget-like request and write resulting file on local disk
import urllib2
request = urllib2.Request('http://ws.resif.fr/fdsnws/dataselect/1/query?network=FR&station=OGDI&channel=HHZ&starttime=2014-01-10T00:00:00&endtime=2014-01-11T00:00:00')
response = urllib2.urlopen(request)
newfile = open('./mydata.mseed', 'w')
newfile.write ( response.read() )

lien complet vers la doc : http://seismology.resif.fr/#CMSConsultPlace:FDSN_WEBSERVICES

En effet comme l'une des demandes du groupe geodata-resif qui va exploiter les logs du téléchargement-prétraitement massif de données est justement d'avoir la taille en octet de ce qui a été téléchargé, l'exemple ci dessus est pile poil ce qu'on veut.

Dans la solution obspy, comme PA l'a tres justement dit ce matin, récupérer cette taille en octet peut se faire grace à l'option filename du get_waveform : https://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_waveforms.html#obspy.clients.fdsn.client.Client.get_waveforms

Quelques notes pour les tests à venir :

  • possibilité de Freeze cache RESIF pour un temps donné
  • option des webservices : non refresh cache de la requête via url (wget ou urllib), mais pas implémenté dans API ObsPy
  • D.W. mentionne la possibilité de configurer wget pour paramétrer le délai et le nombre de retry si échec wget
  • Question : Comment fonctionne le cache obspy en cas de campagnes simultanées ?

-> ToDo pour la prochaine réunion : P.-A.B. + V.L.

3bis)

Liste des codes HTTP rencontrés lors de campagnes CiGri de preprocess avec méthode requêtes wget :

  • HTTP request sent, awaiting response... 200
  • HTTP request sent, awaiting response... 204
  • HTTP request sent, awaiting response... 401
  • HTTP request sent, awaiting response... 500
  • HTTP request sent, awaiting response... 502 Proxy Error
  • HTTP request sent, awaiting response... No data received.
  • HTTP request sent, awaiting response... Read error (Connection reset by peer) in headers.
  • HTTP request sent, awaiting response... 200 OK
  • HTTP request sent, awaiting response... 204 No Content
  • HTTP request sent, awaiting response... 400 Bad Request
  • HTTP request sent, awaiting response... 401 Unauthorized

4)

Pas bon de faire moyenne des Δt en prenant successivement toutes les stations pour référence (car on distribue l'erreur au sein d'un réseau). Ca peut être OK si N très grand 'erreur en 1/N). Plutôt explorer la densité de mesure à chaque pas de temps, parmi les N ref : une branche principale doit ressortir, mettant en évidence un ensemble de référence possible (hypothèse que peu de stations dérivent).

Papier Hable et al : mélange Terre/Mer : pb de distribuer l'erreur sur toutes les stations, stations qui ont des caractéristiques très différentes.

les deux mesures δt positif et δt négatif (avec donc deux cohérences) : pas conclu...

plusieurs Δf pour valider δt ?

dérive dans données synthétiques :

  • créneau : oui
  • linéaire : oui
  • sinusoidale de moyenne nulle : non

le choix de la période de temps pour la réf :

  • J1 uniquement ? -> chaque nouvelle station, sa ref est son J1
  • toute la période ? -> si dérive pendant toute la période, pb : solution est peut être d'itérer sur la mesure du δt (en ajustant progressivement la réf)

choix de la période de référence : dépend du style de déploiement.

Workshop SERBIE Maths Applis: (https://gricad-gitlab.univ-grenoble-alpes.fr/picardch/seismic-monitors-time-delay)

  • 7 stations sur 80 permettent d'avoir un graphe complet pendant 1 an
  • analyze bruit, délai, SNR
  • trouve une structure dans le bruit, pas un bruit blanc

Références de correction de la dérive d'horloge dans les réseaux (systèmes embarqués, ioT) : à compléter ... : non complet, références ou ils modélisent l'horloge

Approche minimisation fonctionnelle + poids sur aretes du graphe

distance normalisée entre les stations [0 1]

cadre analytique avec hyp Δt(k)=0

avantages de l'approche minimisation : fonctionne avec graphe non complet, inclue les cohérence

pyramide de graphe : à l'intérieur d'un graphe incomplet, isoler un graphe complet, déterminer son noeud équivalent (équivalence des graphes), on retrouve un nouveau graphe total qui peut devenir complet.

-> ToDo pour prochaine réunion : Christophe P.

5)

5.1) Argentière

données sur SUMMER dans répertoire RESIF Nodes

mais on laisse un peu en attente la mesure de pseudo-δt Argentière, d'abord finir de traiter données synthétiques + traiter test San Jacinto + traiter test Pyrope

-> ToDo pour prochaine réunion : P.R. explore les données Argentière : bande de fréquence et cie...

5.2) Données synthétiques :

  • nstep = 720
  • SNR=100
  • coupe période en 3, la période du milieu subit une dérive linéaire de 1sec au total : appliquer les dérives sur GF et pas sur δt
  • aussi un créneau
  • ref sur le premier T/3 et ref sur tout T (pb peut etre résolu via itération sur mesure δt)
  • regarder évolution des Δt(t) pour chaque 79 ref, voir la branche principale (représentation en histogramme de évolution des Δt(t,N) de la station qui dérive, à chaque pas de temps (f(t)), et en prenant successivement chaque station pour référence (f(N)).
  • + tester le cas avec plusieurs itérations pour affiner mesure de δt à un pas de temps donné (attention ne pas confondre itération et pas de temps...), si la dérive linéaire se produit sur toute la période de référence. En effet on se repose la question des itérations pour la mesure du δt ? Schönfelder faisait 2 it. Hable en fait autant que nécessaire (jusqu'à 4). Nous on n'a pas trop exploré pour l'instant. Ce sera peut être une solution au cas ou la dérive a lieu pendant toute la période de référence. A explorer donc sur données synthétiques. le δt total est la somme des δt à chaque itération. Par contre à chaque it on aura un nouvel indicateur de cohérence. Il faudrait tout de même que le graphe ne soit pas modifié à chaque itération.

-> ToDo pour prochaine réunion : A.L.

5.3) Pyrope : pour 1ère validation méthodo synthétique

L.S. a mesuré des δt Pyrope

Pour mémoire, on avait déjà cherché à exploiter ses mesures dans cette section : https://ciment.ujf-grenoble.fr/wiki/index.php/Projects/pr-iworms#Calcul_des_.CE.94t_et_Tests_avec_les_profils_PE.2FPF_du_r.C3.A9seau_X7_.28Pyrope.29

mais graphe incomplet : Fonction de Green entre la grille et un profil

L.S. (avec l'aide de Jacques Brives) va rechercher les données pyrope original pour tenter de reconstruire un graphe complet de δt. Voir les stats de disponibilité des données.

-> ToDo pour prochaine réunion : L.S.

5.4) San Jacinto : pour 2ème validation méthodo synthétique

on reprend les correlations 10min de l'EXP003 2.9-5.8Hz (voir wiki San Jacinto : https://ciment.ujf-grenoble.fr/wiki/index.php/Projects/pr-sanjacinto#Production:_S.C3.A9rie_d.27exp.C3.A9riences_00_:_Correlation_Lag_2sec)

Info : V = 450m/s (pour zoom sur onde de surface pour mesure δt) : mais dans un 1er temps on va mesurer δt sur toute la corrélation avec lag +/- 2sec (401 points). En fonction des résultat, on verra dans un 2ème temps pour tenter de séparer les parties causale et acausale.

Pour l'instant on a déjà les corr_10min pendant 1 mois, ainsi que le stack mensuel (ainsi que les 31 stacks journaliers). On propose de faire des stacks horaires avec overlap 10min. On serait donc capable de mesurer des δt et cohérence associée à une fréquence de 10min pendant 1 mois. On va regarder en détail ce qui se passe en milieu de mois (lors du changement de batteries sur les appareils) : y a t il eu des δt associés ?

-> ToDo pour prochaine réunion : A.L.


Ordre du jour :

0) Vacances de Marmitte (mascotte du projet I-Worms):

Marmitte Mariette Soupir
Marmitte.jpg
Mariette.jpg
Soupir.jpg

1) Point d'info prétraitement interp Lanczos (Albanne)

cf plots en section 2 de la réunion du 10/7/2018 : https://ciment.ujf-grenoble.fr/wiki/index.php/Projects/pr-iworms#10.2F07.2F2018_:_R.C3.A9union_14h_en_258_IMAG

2) Point d'info logging prétraitements (Gregory et Jonathan et Catherine et Albanne)

voir réunion 28/6/2018 : https://ciment.ujf-grenoble.fr/wiki/index.php/Projects/pr-iworms#28.2F06.2F2018_:_R.C3.A9union_Logging_pr.C3.A9traitements_10h_245_ISTerre

3) Point d'info code MPI4Py pour prétraitements (P.-A.) -> reporté

3bis) Rapatrie+Preprocess : 3.1) wget/obspy, 3.2) flag quality=None/M (Albanne)

Quality 2018 sept.png
Quality ref.png
QualityM ref.png
Flag Quality None.png
Flag quality M.png

4) Point d'info Summerschool δt -> Δt (Christophe) + Données synthétiques -> vers un papier méthodo (tous)

4.1) Formule analytique graphe complet :

Hypothèse: Δt(N)=0

Solution: pour k=1,...,N-1: Δt(k) = 1/N ( ∑ δtk,i; i=1,N; i≠k + ∑ δti,N; i=1,N-1)

Si on fait toutes les hyp successivement et qu'on moyenne, on obtient la solution :

pour k=1,...,N: Δt(k) = 1/N ( ∑ δtk,i; i=1,N; i≠k ), donc la même solution que dans Hable et al GJI juin 2018 (données RHUM-RUM)

4.2) Plots à synthétiser : exemple ci dessous (pour un SNR donnée, un nb de stations donné, un nb de step dans la réf donné) :

Δti si l'on prend successivement les 79 stations comme référence [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1, SNR=10]
Dt allref gaussianfit SNR 10 SW2 decGF1.png
Dt allref gaussianfit3D SNR 10 SW2 decGF1.png
Distribution et fit gaussien des Δti obtenus en moyennant les N Δti1,i2 provenant du calcul successif en prenant chaque station pour référence [méthode filter pour random noise ; mesure δt avec méthode rapport des phases (CC=0), zoom sur SW avec fenêtre largeur fixe 40sec (SW=2), decGF=1]
Dt R8 gaussianfit meanref SW2 decGF1.png
Dt gaussianfit meanref SW2 decGF1.png

Ces plots peuvent se synthétiser :

SigmaDt.png

5) Quelles données réelles pour test grandeur nature ? SanJacinto, SISMOB Argentière, ... (tous)

5.1) San Jacinto :

On a l'EXP003 qui est dispo (voir https://ciment.ujf-grenoble.fr/wiki/index.php/Projects/pr-sanjacinto#Production:_S.C3.A9rie_d.27exp.C3.A9riences_00_:_Correlation_Lag_2sec)

  • gapfilling : fill remaning gap with zero
  • bandpass filter : 0.5 - 25Hz
  • decimation factor 5 using spline interpolation : -> 100Hz
  • whitening 2.9-5.8Hz with Deltaf = 0.15Hz
  • 1Bit (luke4:/scratch_md1200/lecointre/SAN_JACINTO/FINAL_PREPROCESS_EXP003)
  • compute 10min 100Hz correlations with lag +/- 2sec, for the entire month of data : luke4:/scratch_md1200/lecointre/SAN_JACINTO/CORRELATION/EXP003/SG.2014${JD}.000000.h5
  • compute dailymeans and monthly mean : luke4:/scratch_md1200/lecointre/SAN_JACINTO/CORRELATION/EXP003/DAILYMEAN/SG.DAILYMEAN.2014${JD}.000000.h5 luke4:/scratch_md1200/lecointre/SAN_JACINTO/CORRELATION/EXP003/GLOBALMEAN/SG.GLOBALMEAN.2014.128.158.000000.h5

OK si on moyenne 6 * 10min corr pour mesure δt horaire ?

Plutot que de repartir du hourly bandpass 2-6Hz non 1Bit dispo actuellement pour Chloé  ? (voir https://ciment.ujf-grenoble.fr/wiki/index.php/Projects/pr-sanjacinto#Other_initial_preprocessing) : qui nécessiterait donc 3 étapes supplémentaires :

  • 1-Bit
  • calcul hourly corr
  • compute monthly mean

5.2) SISMOB Argentière :

lecointre@luke:/summer$ sg pr-resolve -c "ls resolve/Argentiere/DATA/"
GPR	  KEPHREN-SISMO	RESIF_Nodes_01	      SISMO
GPS_RTK  RESIF_Nodes	RESIF_Nodes_20180925


enregistrement continu de 100 capteurs 3C pendant un peu plus de 30 jours.

dérive en distance relative de la nappe de capteurs au cours du temps.

quand des capteurs se mettent à dysfonctionner ou que les corrélations ne marchent plus bien : apport de techniques de machine learning, sparcity/compacity ?

rencontre avec Olivier Michel (GIPSA-Lab) prévue le 5/10

[edit] 22-10-2018 : Reunion preprocess MPI4Py

Participants : P.-A. Bouttier, A. Lecointre

Durée : 30min

Compte-Rendu :

[edit] 30-11-2018 : Reunion 10h-12h IMAG 248

Participants : P.A. Bouttier (30'), V. Louvet, P. Roux, J. Schaeffer, F. Thollard, A. Lecointre

Durée : 1h45

Compte-Rendu :

0) Déploiement MPI4Py du get_input + preprocess

utilise la liste de 237 stations (traces) suivantes : File:StationsByCodes.txt File:List rsync.txt obtenues avec les params suivants :

JJJ=2017179
listnet="FR,RD,RA,ZH,Z3,G"
liststa="*"
listloc="*"
listchan="HH*"
minlat=43
maxlat=46
minlon=3
maxlon=8

Prêt à être déployé sur luke

Bémol : bug OGS3 2017179 interpolation avec new starttime: Si l'on active new starttime, et si l'on est dans l'env P.-A. (nix, python3, mpi4py, ...), l'interpolation de la 4ème et dernière trace (itr=3) ne passe pas :

Exemple dans le contexte ou ca ne bugge pas (python2.7, environement module, séquentiel)

('itr:', 0)
4 Trace(s) in Stream:
FR.OGS3.00.HHN | 2017-06-27T23:58:47.955000Z - 2017-06-28T00:50:01.700000Z | 200.0 Hz, 614750 samples
FR.OGS3.00.HHN | 2017-06-28T00:50:20.245000Z - 2017-06-28T11:06:38.380000Z | 200.0 Hz, 7395628 samples
FR.OGS3.00.HHN | 2017-06-28T11:07:07.225000Z - 2017-06-28T11:27:34.980000Z | 200.0 Hz, 245552 samples
FR.OGS3.00.HHN | 2017-06-28T11:27:37.045000Z - 2017-06-29T00:01:17.020000Z | 200.0 Hz, 9043996 samples
2017-06-27T23:58:47.960000Z
[1026 1023 1019 ..., 1016 1018 1021]
('itr:', 1)
[...]
('itr:', 3)
4 Trace(s) in Stream:
FR.OGS3.00.HHN | 2017-06-27T23:58:47.960000Z - 2017-06-28T00:50:01.690000Z | 100.0 Hz, 307374 samples
FR.OGS3.00.HHN | 2017-06-28T00:50:20.250000Z - 2017-06-28T11:06:38.370000Z | 100.0 Hz, 3697813 samples
FR.OGS3.00.HHN | 2017-06-28T11:07:07.230000Z - 2017-06-28T11:27:34.980000Z | 100.0 Hz, 122776 samples
FR.OGS3.00.HHN | 2017-06-28T11:27:37.045000Z - 2017-06-29T00:01:17.020000Z | 200.0 Hz, 9043996 samples
2017-06-28T11:27:37.050000Z
[ 978  976  972 ..., 1051 1048 1045]
---> endroit du bug dans l'autre contexte d'exécution <---
4 Trace(s) in Stream:
FR.OGS3.00.HHN | 2017-06-27T23:58:47.960000Z - 2017-06-28T00:50:01.690000Z | 100.0 Hz, 307374 samples
FR.OGS3.00.HHN | 2017-06-28T00:50:20.250000Z - 2017-06-28T11:06:38.370000Z | 100.0 Hz, 3697813 samples
FR.OGS3.00.HHN | 2017-06-28T11:07:07.230000Z - 2017-06-28T11:27:34.980000Z | 100.0 Hz, 122776 samples
FR.OGS3.00.HHN | 2017-06-28T11:27:37.050000Z - 2017-06-29T00:01:17.010000Z | 100.0 Hz, 4521997 samples
[  975.99390635   970.99872303   971.99992675 ...,  1047.00000655
 1044.00006954  1050.9973736 ]

Exemple dans le contexte ou ca bug :

[...]
3
4 Trace(s) in Stream:
FR.OGS3.00.HHN | 2017-06-27T23:58:47.960000Z - 2017-06-28T00:50:01.700000Z | 100.0 Hz, 307375 samples
FR.OGS3.00.HHN | 2017-06-28T00:50:20.250000Z - 2017-06-28T11:06:38.380000Z | 100.0 Hz, 3697814 samples
FR.OGS3.00.HHN | 2017-06-28T11:07:07.230000Z - 2017-06-28T11:27:34.980000Z | 100.0 Hz, 122776 samples
FR.OGS3.00.HHN | 2017-06-28T11:27:37.045000Z - 2017-06-29T00:01:17.020000Z | 200.0 Hz, 9043996 samples
2017-06-28T11:27:37.050000Z
[ 978  976  972 ... 1051 1048 1045]
---> endroit du bug <---

On a cette erreur :

Traceback (most recent call last):
 File "./preprocess_refact.py", line 272, in <module>
   s = interpol(s, int(syr), int(sjd))
 File "./preprocess_refact.py", line 104, in interpol
   s[itr].interpolate(sampling_rate=ff,starttime=nts,method='lanczos',a=5)
 File "<decorator-gen-27>", line 2, in interpolate
 File "/home/lecointre/iworms/iworms_preprocess_src/version_nix/virtualenv/obspy/lib/python3.6/site-packages/obspy/core/util/decorator.py", line 246, in skip_if_no_data
   return func(*args, **kwargs)
 File "<decorator-gen-26>", line 2, in interpolate
 File "/home/lecointre/iworms/iworms_preprocess_src/version_nix/virtualenv/obspy/lib/python3.6/site-packages/obspy/core/util/decorator.py", line 235, in raise_if_masked
   return func(*args, **kwargs)
 File "<decorator-gen-25>", line 2, in interpolate
 File "/home/lecointre/iworms/iworms_preprocess_src/version_nix/virtualenv/obspy/lib/python3.6/site-packages/obspy/core/trace.py", line 232, in _add_processing_info
   result = func(*args, **kwargs)
 File "/home/lecointre/iworms/iworms_preprocess_src/version_nix/virtualenv/obspy/lib/python3.6/site-packages/obspy/core/trace.py", line 2412, in interpolate
   starttime, dt, npts, type=method, *args, **kwargs))
 File "/home/lecointre/iworms/iworms_preprocess_src/version_nix/virtualenv/obspy/lib/python3.6/site-packages/obspy/signal/interpolation.py", line 291, in lanczos_interpolation
   _validate_parameters(data, old_start, old_dt, new_start, new_dt, new_npts)
 File "/home/lecointre/iworms/iworms_preprocess_src/version_nix/virtualenv/obspy/lib/python3.6/site-packages/obspy/signal/interpolation.py", line 39, in _validate_parameters
   raise ValueError("The new array must be fully contained in the old "
ValueError: The new array must be fully contained in the old array. No extrapolation can be performed.
Traceback (most recent call last):
 File "h5py/_objects.pyx", line 193, in h5py._objects.ObjectID.__dealloc__
RuntimeError: Can't decrement id ref count (can't close file, there are objects still open)
Exception ignored in: 'h5py._objects.ObjectID.__dealloc__'
Traceback (most recent call last):
 File "h5py/_objects.pyx", line 193, in h5py._objects.ObjectID.__dealloc__
RuntimeError: Can't decrement id ref count (can't close file, there are objects still open)

C'est intéressant de constater que même pour les 2 premières traces, il y a une petite différence (à la fin ca coupe pas tout à fait au même point).

Pour avancer, on propose de zapper l'interpolation, ou de traiter l'exception.

On peut ajouter des dates de l'été 2016 à l'été 2017 : de 2016183 à 2017181 (366 jours).

On fixe une session de travail le 13/12 à 9h à ISTerre en présence des ingés RESIF pour déployer les tests de montée en charge RESIF progressive. Laisser actif cache RESIF applicatif. Cache luke aussi. Réservation noeuds dahu pour tests -> P.-A. s'en occupe ?

Ensuite la numérotation des points reprend la ToDo list de réunion précédente du 01/10

5.2) Données synthétiques :

  • nstep = 720
  • SNR=100
  • coupe période en 3, la période du milieu subit une dérive linéaire de 1sec au total : appliquer les dérives sur GF et pas sur δt
  • aussi un créneau
  • ref sur le premier T/3 et ref sur tout T (pb peut etre résolu via itération sur mesure δt : en fait NON)
  • regarder évolution des Δt(t) pour chaque 79 ref, voir la branche principale (représentation en histogramme de évolution des Δt(t,N) de la station qui dérive, à chaque pas de temps (f(t)), et en prenant successivement chaque station pour référence (f(N)).
  • + tester le cas avec plusieurs itérations pour affiner mesure de δt à un pas de temps donné (attention ne pas confondre itération et pas de temps...), si la dérive linéaire se produit sur toute la période de référence. En effet on se repose la question des itérations pour la mesure du δt ? Schönfelder faisait 2 it. Hable en fait autant que nécessaire (jusqu'à 4). Nous on n'a pas trop exploré pour l'instant. Ce sera peut être une solution au cas ou la dérive a lieu pendant toute la période de référence en fait NON. Le δt total est la somme des δt à chaque itération. Par contre à chaque it on aura un nouvel indicateur de cohérence. Il faudrait tout de même que le graphe ne soit pas modifié à chaque itération ?
Station 1 dérive (uniquement S1)
niter=1 Δt niter=5 Δt niter=5 δt
Tref=[0:T/3], creneau 1sec
Synt Dt 379483.png
Synt Dt 379535.png
Synt dtbypair 379535.png
Tref=[0:T], creneau 1sec
Synt Dt 379487.png
Synt Dt 379537.png
Synt dtbypair 379537.png
Tref=[0:T/3], lineaire 1sec
Synt Dt 379482.png
Synt Dt 379533.png
Synt dtbypair 379533.png
Tref=[0:T], lineaire 1sec
Synt Dt 379486.png
Synt Dt 379534.png
Synt dtbypair 379534.png

Action : rajouter des plots avec :

  • Dt station1 en moyennant sur toutes les stations de référence --> OK Ajout courbe tiretés rouge
  • Dt station2 (qui ne dérive pas) --> OK ajout colonne 1 tableau ci dessous, uniquement pour niter=1
  • appliquer une dérive à 2 stations au lieu d'une seule --> OK ajout colonnes 2-3 tableau ci dessous, uniquement pour niter=1
  • test aussi avec SNR=1 --> OK ajout colonne 4 tableau ci dessous, uniquement pour niter=1 et seule S1 dérive
niter=1 Δt
SNR=100, Seule S1 dérive (plot S2). SNR=100, S1 et S2 dérivent (plot S1). SNR=100, S1 et S2 dérivent (plot S2). SNR=1, Seule S1 dérive (plot S1).
Tref=[0:T/3], creneau 1sec
Synt Dt 379483 2.png
Synt Dt 451846.png
Synt Dt 451846 2.png
Synt Dt 451863.png
Tref=[0:T], creneau 1sec
Synt Dt 379487 2.png
Synt Dt 451856.png
Synt Dt 451856 2.png
Synt Dt 451865.png
Tref=[0:T/3], lineaire 1sec
Synt Dt 379482 2.png
Synt Dt 451855.png
Synt Dt 451855 2.png
Synt Dt 451864.png
Tref=[0:T], lineaire 1sec
Synt Dt 379486 2.png
Synt Dt 451857.png
Synt Dt 451857 2.png
Synt Dt 451866.png

5.3) Pyrope Maupasacq : pour 1ère validation méthodo synthétique

Maupasacq: 184 pas de temps, 1 pas de temps = 1 jour

Correlations sur 3601 points : f=5Hz, maxlag=360sec : 360*5*2+1 = 3601 points.

Quelle fréq centrale de spectre ? Au hasard 4.35 Hz (comme SJ)...

441 stations -> 97461 paires (=441*442/2 donc avec autocorr, ex les indices 0 et 5 sont des autocorr : dt ~ 1E-15 et corr ~ 0.90-0.92). Attention les paires ne sont pas ordonnées.

dt coh
Dt maupasacq.png
Coh maupasacq.png

En fait le graphe n'est pas complet : Nan Values : soit std(GF_timestep)=0, soit std(GFref)=0 :

lecointre@luke:~/iworms/iworms_measure_dt_src/MAUPASACQ$ grep GFref OAR.measure_dt.8459270.stdout|wc -l
443
lecointre@luke:~/iworms/iworms_measure_dt_src/MAUPASACQ$ grep GFref OAR.measure_dt.8459270.stdout|head -5
(1818, u'XD.N0102.00_xD.MBB17.00_ZZ', 0.01176309585571289, 'np.std(GFref)=0')
(2137, u'XD.N0103.00_xD.MBB17.00_ZZ', 0.0016770362854003906, 'np.std(GFref)=0')
(2460, u'XD.N0105.00_xD.MBB17.00_ZZ', 0.008351802825927734, 'np.std(GFref)=0')
(2830, u'XD.N0106.00_xD.MBB17.00_ZZ', 0.013605117797851562, 'np.std(GFref)=0')
(3126, u'XD.N0108.00_xD.MBB17.00_ZZ', 0.0016019344329833984, 'np.std(GFref)=0')

Action:

  • voir avec Laurent pourquoi graphe incomplet Maupasacq (GF constante ?) : Il y a effectivement des GF=0 pour certaines paires et certains pas de temps. D'après Pierre Boué, leur méthode de calcul des corrélations force GF=0 dès lors que la corrélation est suspecte. On laisse tomber ces données incomplètes...
  • insister pour avoir plutot Pyrope graphe complet : Non, on laisse tomber ces données incomplètes.
  • regarder la distribution des dt et coh associée, voir comment ca fitte avec les tests synthétiques : On mesure des δt max = 0.1149 : ci dessous :
distribution des dt et coh
Dt distribution Maupasacq.png

F=5Hz donc la cohérence associée est en fait toujours celle pour $δt = 0 ?

La cohérence associée est faible ? rien à voir avec δt mesuré sur données synthétiques ???

5.4) San Jacinto : pour 2ème validation méthodo synthétique

mesure de dt et coh, toutes les 10min pendant 2 jours (JD=140 et JD=154), par rapport à une référence mensuelle. Fcentrale = 4.35 Hz (EXP003)

dt coh
JD=140
Dt SJ 140.png
Coh SJ 140.png
JD=154
Dt SJ 154.png
Coh SJ 154.png

Je confirme qu'il y a une inversion entre les axes abscisse et ordonnée sur ces figures SJ (Par contre les fig maupasacq au dessus sont OK)

Actions:

  • regarder la distribution des dt et coh associée, voir comment ca fitte avec les tests synthétiques
  • tester l'inversion dt -> Dt (graphe complet, avec une station pour ref, avec toutes les stations prises successivement pour ref)
  • regarder la distribution des Dt, voir comment ca fitte avec les tests synthétiques

[edit] 06-12-2018 : Reunion 10h ISTerre : discussion logging

Participants : G. Arneodo, A. Lecointre

Durée : 30min

Compte-Rendu :

On se met d'accord sur un format de log de sortie pour le preprocess, qui soit exploitable avec ELK (Elasticsearch - Logstash - Kibana).

Albanne instrumente le code preprocess_refact.py : chaque processus MPI écrit dans un fichier ascii dédié, puis, après que mpirun ait rendu la main, le batch submit script concatène ces fichiers ascii en un seul log file.

Exemple de log file File:Log.txt final (exécuté sur 10 coeurs luke4) : (stderr vide et stdout: File:OAR.logging.9543481.stdout.txt)

date hostname userid net sta loc chan wanted_starttime wanted_endtime nb_downloaded_samples total_elapsed_time:rapatrie+preprocess+writeintoh5file

lecointre@luke:~/iworms/iworms_preprocess_src/version_nix/LOG$ head -15 log.txt
2018-12-07+15:40:50Z luke4 lecointrea FR ARBF 00 HHE  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8191567 8.428675
2018-12-07+15:40:59Z luke4 lecointrea FR BLAF 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 6334126 6.668895
2018-12-07+15:41:06Z luke4 lecointrea FR CFF 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8656993 41.807562
2018-12-07+15:41:47Z luke4 lecointrea FR ENAUX 00 HHE  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8652766 41.352341
2018-12-07+15:42:29Z luke4 lecointrea FR GRN 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 4227324 25.149063
2018-12-07+15:42:54Z luke4 lecointrea FR MLYF 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 1094433 9.017803
2018-12-07+15:43:03Z luke4 lecointrea FR OGCB 00 HHE  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8657088 41.515094
2018-12-07+15:43:44Z luke4 lecointrea FR OGMO 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8656250 41.520346
2018-12-07+15:44:26Z luke4 lecointrea FR OGS2 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 17284572 42.347195
2018-12-07+15:45:08Z luke4 lecointrea FR OGVG 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8657312 41.578747
2018-12-07+15:45:50Z luke4 lecointrea FR RSL 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 6890368 6.564258
2018-12-07+15:45:56Z luke4 lecointrea FR RUSF 05 HHE  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 0 4.457264
2018-12-07+15:46:01Z luke4 lecointrea FR SAOF 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 6301048 6.500160
2018-12-07+15:46:07Z luke4 lecointrea FR SURF 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 6567600 6.309196
2018-12-07+15:46:14Z luke4 lecointrea FR VAL4 00 HH1  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 0 4.472590
lecointre@luke:~/iworms/iworms_preprocess_src/version_nix/LOG$ tail -5 log.txt 
2018-12-07+15:47:01Z luke4 lecointrea Z3 A192A 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8655978 41.474784
2018-12-07+15:47:42Z luke4 lecointrea Z3 A196A 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8654746 41.467375
2018-12-07+15:48:24Z luke4 lecointrea Z3 A202A 00 HHE  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8653602 41.706700
2018-12-07+15:49:05Z luke4 lecointrea Z3 A206A 00 HHN  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8655920 41.483002
2018-12-07+15:49:47Z luke4 lecointrea Z3 A217A 00 HHZ  2017-06-27+23:59:00Z 2017-06-29+00:01:00Z 8648736 41.466532
lecointre@luke:~/iworms/iworms_preprocess_src/version_nix/LOG$ wc -l log.txt 
236 log.txt

Les traces pour lesquelles il n'y avait pas de données apparaissent avec un 0 dans le champ nb_downloaded_samples:

lecointre@luke:~/iworms/iworms_preprocess_src/version_nix$ grep "Z 0 " LOG/log.txt |wc -l
18
lecointre@luke:~/iworms/iworms_preprocess_src/version_nix$ grep "This trace won' be processed" OAR.logging.9543481.stdout|wc -l
18

[edit] 13-12-2018 : Session de travail 9h-12h ISTerre : deploiement preprocess

Participants : V. Louvet, P.-A. Bouttier, G. Arneodo, A. Lecointre (+ D. Woliniec, J. Schaeffer, R. Bouazzouz, C. Pequegnat)

Durée : 3h

Compte-Rendu :

code déployé sur 2-3 noeuds froggy (16cores/node) et/ou dahu (32cores/node).

sur Froggy ca freeze au delà de 4 nodes : je ne me souviens plus si on a compris pourquoi ?

sur Dahu : serveur RESIF rejette nos requêtes à partir de 3 nodes (=96 MPI processes). Gregory confirme la présence d'une limitation à 80 requêtes simultanées de leur côté.

Gregory augmente sa limite à 120.

On parvient à tourner sur 3 (voire 4) noeuds dahu.

Mais on atteint les limites CPU de la VM qui supporte le serveur RESIF (donc cette limite à 80 était finalement bien dimensionnée pour les capacités de la machine).

Il reste des questions en suspens :

1) comprendre traiter proprement les erreurs d'interpolation de ce type (qui ne se produisent pas que pour OGS3 HHN) :

Traceback (most recent call last):
 File "./preprocess_refact.py", line 265, in <module>
   s = interpol(s, int(syr), int(sjd))
 File "./preprocess_refact.py", line 99, in interpol
   s[itr].interpolate(sampling_rate=ff,starttime=nts,method='nearest')
 File "<decorator-gen-27>", line 2, in interpolate
 File "/Users/pierre-antoinebouttier/.miniconda3/envs/obspy/lib/python3.6/site-packages/obspy/core/util/decorator.py", line 246, in skip_if_no_data
   return func(*args, **kwargs)
 File "<decorator-gen-26>", line 2, in interpolate
 File "/Users/pierre-antoinebouttier/.miniconda3/envs/obspy/lib/python3.6/site-packages/obspy/core/util/decorator.py", line 235, in raise_if_masked
   return func(*args, **kwargs)
 File "<decorator-gen-25>", line 2, in interpolate
 File "/Users/pierre-antoinebouttier/.miniconda3/envs/obspy/lib/python3.6/site-packages/obspy/core/trace.py", line 232, in _add_processing_info
   result = func(*args, **kwargs)
 File "/Users/pierre-antoinebouttier/.miniconda3/envs/obspy/lib/python3.6/site-packages/obspy/core/trace.py", line 2412, in interpolate
   starttime, dt, npts, type=method, *args, **kwargs))
 File "/Users/pierre-antoinebouttier/.miniconda3/envs/obspy/lib/python3.6/site-packages/obspy/signal/interpolation.py", line 72, in interpolate_1d
   new_start, new_dt, new_npts)
 File "/Users/pierre-antoinebouttier/.miniconda3/envs/obspy/lib/python3.6/site-packages/obspy/signal/interpolation.py", line 39, in _validate_parameters
   raise ValueError("The new array must be fully contained in the old "
ValueError: The new array must be fully contained in the old array. No extrapolation can be performed.
Traceback (most recent call last):
 File "h5py/_objects.pyx", line 195, in h5py._objects.ObjectID.__dealloc__
RuntimeError: Can't decrement id ref count (Can't close file, there are objects still open)
Exception ignored in: 'h5py._objects.ObjectID.__dealloc__'
Traceback (most recent call last):
 File "h5py/_objects.pyx", line 195, in h5py._objects.ObjectID.__dealloc__
RuntimeError: Can't decrement id ref count (Can't close file, there are objects still open)

2) comprendre ce que signifient ces exceptions : "No FDSN service ..." a t on affaire au même type de problème qu'ici (déjà rencontré lors du WGET rapatrie + traitement sur grille cigri)  : doit on faire des retry ? combien ? avec quel délai ?

/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: (null): Warning: Data integrity check for Steim-1 failed, last_data=-2292332, xn=-2325845
warnings.warn(*_i)
/applis/ciment/v2/stow/x86_64/gcc_4.6.2/python_2.7.12/lib/python2.7/site-packages/obspy/io/mseed/core.py:423: InternalMSEEDReadingWarning: readMSEEDBuffer(): Last record only has 63 byte(s) which is not enough to constitute a full SEED record. Corrupt data? Record will be skipped.
warnings.warn(*_i)

La grosse difficulté est que ces erreurs sont difficiles à reproduire.

[edit] 21-01-2019 : Reunion 9h30 IMAG 258

Participants :

Durée :

Compte-Rendu :

[edit] 04 et 05-02-2019 : Workshop du projet, Chamonix

Participants : PAB, PR, CP, VL, AL

Durée : 2jours

Compte-Rendu :

[edit] 21-01-2019 : Reunion 10h00 IMAG 258

Participants : PAB, VL, CP, PR, AL

Durée : 1h

Compte-Rendu :

1. Papier :

  • à rédiger partie formelle méthodo dt -> Dt dans le cas du graphe complet + pistes pour graphes incomplet, mais pas de mise en pratique
  • à relire partie challenge et métriques du calcul : chiffres.
  • à tester jupyter notebook : synthetic.ipynb

2. San Jacinto

conclusions du scatter plot Dt,Coh à rédiger dans papier:

  • réseau bien synchronisé
  • Dt -> Dtw : on resserre autour de zéro : pondération par cohérence permet d'affaiblir les fortes amplitudes de dt.
  • valeurs extrêmes Dt = 0.05sec -> freq centrale = 4.35Hz -> 0.23sec de période -> 0.05/0.23 = 22% de période

movie Dt:

  • cohérence spatiale des Dt, alors que rien dans la mesure du dt puis calcul du Dt ne fait intervenir un lissage spatial

comparaison météo:

  • dans un 1er temps : FFT météo et FFT(Dt) -> corr sur spectre
  • voir si différentes fréquences se relient à différents proxy météo : 24h sur T° , plus haute fréq sur vent, ... ?
  • sur station 846 (fort coh et faible Dt)
  • peut etre effacer préalablement les 2 jours avec faible coh (LUNDI 139 et MARDI 147)

comparaison carte de vitesse à 200m:

  • comparer )à std temporelle de Dt
  • comparer à moy temporelle de coh (sans les 2 jours problématiques)

moments : espérance, variance, ...

piste : wavelet transforms (spatial 2D en x,y) sur V_200m(x,y) et Dt(x,y,t) -> à corréler entre elles

plots dynamiques "octave" par box (dx,dy) -> qqch à relier avec humidité en milieu de période

San Jacinto : 2 refs à pointer dans le papier :

  • papier Yehuda pour description manip (déploiement capteurs réseua dense, 1 mois, faille, ...)
  • papier Roux GJI pour calcul des corrélations (utilisées pour mesure dt)

[edit] 04-06-2019 : Reunion 14h00 IMAG 258

Participants : PAB, CP, PR, AL

Durée : 1h

Compte-Rendu :

Papier :

  • Ajout ref S. Hable dans intro
  • spectre Dtw a chaque station -> saturation colorbar 1e-3 -> 1e-5 (facteur 100 only)

-> ajout moyenne (spectre moyen) des 1108 stations

-> corréler le spectre de chaque station avec les forcage météo : représenter la station qui fitte le mieux -> illustre le lien entre forcage environnementaux (meteo) et valeur de Dt pour la station

-> separer les spectres des différents indicateurs météo

  • 1d et ses sous harmoniques (12h, 6h) -> temperature ?
  • 8h -> durée de la nuit ?
  • et le 9h ???

---

  • Mordret : coupe vitesse à 200m : voir les autres profondeurs. si on va plus en surface ?
  • retirer les 2 jours pathologiques (en coh : donc représenter moyenne_sur_1108_station(coh) = f(t)) : 139 et 147 : MAJ figures spatiale mean_temporelle(Coh) et std_temporelle(Dtw) : superposer le contour de la faille
  • lineariser spatialement mean_temp(coh)[:] et std_temp(Dtw)[:] : vecteurs de 1108 éléments : que donne la corr ? a t on bien pour chaque station : std(Dtw) augmente quand mean(Coh) diminue ?
  • plot spatial std(Dtw) : forcer encore saturation : voir aussi les extremes hors saturation : si c'est bien 5msec. pour 4.35Hz -> 5msec / (1/4.35) = 2.175% -> a rlier aux valeurs de mesures de doublets.

2 corréaltions :

  • avec la structure : les plots spatiaux
  • avec les forcages externes : les plots spectre meteo

[edit] Mercredi 25-09-2019 : Reunion 10h00 IMAG 258

Participants : PAB, CP, PR, VL, AL

Durée :

Ordre du jour :

  • question : matrice d'adjacence : zeros sur la diagonale : donc pas d'aretes qui partent et arrivent au meme noeud. or j'ai inclus ces 1 sur la diagonale de coh_2D pour le calcul de Coh et AbsCoh : peut etre devrait on les masker ?
  • question : meme question pour la normalisation du Dtw
  • question : et meme pour le calcul de Dt : est ce que cela signifie que l'on devrait diviser par N-1 au lieu de N ? Mais alors la manip avec le pivot de Gauss n'est plus OK...
  • question : pour ce calcul Dtw : comment justifier les valeurs absolues autour du coh. Cas d'une mesure éventuelle de dt associée à coh=-1 : on lui donnerait le poids de coh=1 ???
  • partie computing challenges : il me manque des chiffres : rejouer un bench sur 2-3 noeuds DAHU (max 100 MPI proc) : install nix h5py-mpi ?
  • partie application san jacinto structure spatiale : trouver la bonne représentation.

voir plots 2D lon,lat : https://gricad-gitlab.univ-grenoble-alpes.fr/lecoinal/iworms/wikis/sanjacinto5

Compte-Rendu :