Light Pollution

La mesure de la pollution lumineuse est une mesure clé pour estimer son impact sur notre biosphère. Ces mesures combinent la photométrie à la spectrométrie, pour déterminer exactement l’influence d’un luminaire sur l’environnement.

Pour connaître la théorie et comprendre le sujet, voici quelques sites :

http://www.ascen.be/
https://www.darksky.org/
https://fr.wikipedia.org/wiki/Pollution_lumineuse
https://www.globeatnight.org/light-pollution.php
http://lightpollution.beastro.be/doku.php

Quels sont les objets mesurés et les capteurs adaptés ?

Mesurer la luminance du ciel nocturne

Elle peut généralement se mesurer au niveau amateur depuis le sol via trois outils :

  • Un “luminomètre”, un appareil généralement assez couteux pour disposer d’une mesure suffisamment sensible.
  • Un sky quality meter (SQM) (ou instrument de mesure de la qualité du ciel) utilisant des unités en “magnitudes par seconde d’arc carré” sur des longueurs d’ondes “visibles”. Accessible financièrement, il existe diverses versions, des “compatibles” et même des “DIY”
  • Des photos du ciel… A partir du moment où on calibre son APN pour le faire.

Exemple 1 : Traitement des données cumulatives issues d’un SQM/LU/DL Unihedron

Ce type de capteur est un “enregistreur” de données automatique, parfaitement indépendant (6 piles = plusieurs mois de données, selon la fréquence). Mais on peut récupérer celles-ci de temps en temps, ce qui mène à disposer de plusieurs fichiers cumulatifs couvrant des périodes identiques.

On doit donc trier, supprimer les doubles et rapidement obtenir des données de tendances. Je vous propose ici deux manières de réaliser cela via programmation : du “pur” Python et Python + Pandas

Le fichier “.dat” fournit par le logiciel contient des lignes utiles pour gérer les informations de localisation (préfixées par “#”), mais ne sont pas utiles ensuite. Je les supprime automatiquement au sauvetage dans la bonne directory de stockage par localisation/date. Ensuite, les données utiles sont stockées en philosophie “comma separated” (mais avec “;”).
Un exemple :

2018-04-08T18:30:10.000;2018-04-08T20:30:10.000;18.0;4.91;9.68
=> UTC; LT; Temp(°C); Volt; Mag (en mag/arcsec²) 

Pur Python : Lire, ordonner, calibrer, supprimer les données inutiles (<=10 mag/arcsec²), enlever les doubles, sauver…

# -*- coding: utf-8 -*-
"""
Created on Sun May 27 09:25:40 2018
@author: TTF
First mandatory step : all SQM csv files are concatenated into a single one, 
by a simple "type *.dat > SQM_Global.csv" command and previously "cleaned" for useless lines at capture time
"""
import os, sys,shutil
import time
#start operation
start=time.time()
#First, sort file
with open("F:\\data\\Logs_SQM\\SQM_Global.csv","r") as f:
    sorted_file = sorted(f)
with open("F:\\data\\Logs_SQM\\SQM_Global_Sort.csv","w") as f:
    f.writelines(sorted_file)
#Now examine each lines  
SQMin = open("F:\\data\\Logs_SQM\\SQM_Global_Sort.csv", "r",encoding="utf-8")
SQMout = open("F:\\data\\Logs_SQM\\PI_SQM_Global.csv", "w",encoding="utf-8")
Previous = ""
Current = SQMin.readline() # read a first record
Calib = 0.11
Count_Lu = 0
Count_Double = 0
Count_Valid = 0
Count_Invalid = 0
Magval = 0.0

while Current != "": # end of file 
    Count_Lu = Count_Lu + 1
    if Current != Previous :
        UTC, LC, Temp, Volt, Mag = Current.split(';')
        Magval = float(Mag)
        if (Magval + Calib) >= 10:
            rec_out = LC + "," + Temp + "," + Volt + "," + str(Magval + Calib) + "\n"        
            SQMout.write(rec_out)
            Count_Valid = Count_Valid + 1
        else:
            Count_Invalid = Count_Invalid + 1
        Previous = Current
        Current = SQMin.readline()
    else:
        Count_Double = Count_Double + 1
        Current = SQMin.readline()    
    
SQMin.close()
SQMout.close()
print("Nombre mesures lues      : "+str(Count_Lu))
print("Nombre mesures en double : "+str(Count_Double))
print("Nombre mesures invalides : "+str(Count_Invalid))
print("Nombre mesures Valides   : "+str(Count_Valid))
end=time.time()
print("Start Time               : "+str(start))
print("End Time                 : "+str(end))
print("Duration                 : "+str(end-start)+" sec")

=> Pour trier et traiter 21625 lignes = 0,05 sec…

Python + Pandas : processus identique…

# -*- coding: utf-8 -*-
"""
Created on Sat May 26 07:29:05 2018
@author: TTF
First mandatory step : all SQM csv files are concatenated into a single one, 
by a simple "type *.dat > SQM_Global.csv" command and previously "cleaned" for useless lines at capture time
"""
from pandas import *
import time
#start operation
start=time.time()
#read all data
sqm_df = read_csv("F:\\data\\Logs_SQM\\SQM_Global.csv",
                  delimiter=';',
                  names=['UTC','LC','Temp','Volt','Mag'],
                  parse_dates={'datetime': ['LC']})
#sort by date and time 
df = sqm_df.sort_values(by='datetime')
#add calibration
sqm_df['Calib'] = 0.11
sqm_df['Mag'] = sqm_df['Mag'] + sqm_df['Calib']
sqmtoclean_df = sqm_df.drop(columns='Calib')
#suppress duplicate records
sqmclean_df = sqmtoclean_df.drop_duplicates(keep='first')
#suppress useless records
sqmdata_df = sqmclean_df[sqmclean_df["Mag"] >= 10]
#save epured data 
sqmout = sqmdata_df.loc[:,['datetime','Temp','Volt','Mag']]
sqmout.to_csv("F:\\data\\Logs_SQM\\Pandas_SQM_Global.csv",
              index=False)
end=time.time()
print("Start Time               : "+str(start))
print("End Time                 : "+str(end))
print("Duration                 : "+str(end-start)+" sec")

Quelques remarques :
– calibration : dans les deux cas, on “ajoute” 0.11 mag/arcsec² à chaque mesure, une différence due à la transparence du plastique qui protège le capteur dans sa protection “all weather” (à déterminer selon le cas…)
– syntaxe et possibilités… Le but n’est pas de faire un cours de Pandas, mais cette librairie est dédiée “Data analysis” et offre des capacités d’analyse statistique étendues (Corrélation, déviation, covariance, Kurtosis, Asymétrie, etc…) en quelques commandes. Par exemple :

from pandas import *
#read all data
sqm_df = read_csv("F:\\data\\Logs_SQM\\Pandas_SQM_Global.csv")
plt.figure(); sqm_df.plot()
En 1 ligne de code…

Ou on découvre qu’il a fait jusque -15°C … 🙂
Et si c’est le cas, il y a-t-il une relation entre température et qualité du ciel ?
Q : s’il fait froid, le ciel est-t-il plus dégagé, donc plus noir ?

from pandas import *
#read all data
sqm_df = read_csv("F:\\data\\Logs_SQM\\Pandas_SQM_Global.csv")
#only keep SQM over 17 Mag/arcsec²
sqmdata_df = sqm_df[sqm_df["Mag"] >= 17]
#analyse data by plotting them
sqmdata_df.plot(kind="scatter",x="Temp", y="Mag")

R : Apparemment, de 0 à -8°, c’est vrai.. Mais a une température plus basse (<10°) correspond une magnitude proche de 19, ce qui serait plutôt l’inverse…!

Exemple 2 : “Automatic archive” de données SQM

Tout le monde n’a pas un serveur avec une DB sous la main… Et d’ailleurs, quand on veut garder des données longtemps, autant avoir des copies “fichiers” quelque part…

Mais d’un autre côté, quand on revient d’une prise de mesure, ou on met à jour (trop) rapidement ses données avec la capture du jour (entre deux autres activités) : “shit can happen” 🙂 Et c’est à ce moment-là que l’on s’aperçoit que l’on a mélange/effacé/mal stocké LE fichier que vous avez justement besoin maintenant ! (Murphy’s Law application)

Mesure : les SQM de Unihedron peuvent être à lecture immédiate (/LE/LU) ou différée (/DL). Dans ce dernier cas, le boitier fonctionne seul sur batterie, les mesures sont prises à une période fixe (ex : 30 min) et les données sont stockées, à localement (jusque 32K mesures, donc à 30 Min/mesure = 2 ans de données !). Si on efface la mémoire locale, mieux vaut ne pas perdre ses fichiers ! Mais.. La connexion pour les récupérer est lente (“série” sur USB) et si on laisse les mesures dans la mémoire, elles se cumulent à chaque sauvetage => exploitation immédiate difficile.

Donc, le programme suivant vise à minimiser les erreurs humaines…

Son but est d’archiver automatiquement les fichiers de données (sans trop faire attention au contenu cumulé ou pas), avec des “versions” de l’archive (en cas de problèmes : on peut repartir avec la version précédente).

Seule obligation : standardiser les fichiers de capture selon le format suivant :
<date>_<heure>_<location>(_<calibration>).dat
<date> : format yyyymmdd
<heure> : format hhmmss
<location> : une chaîne de caractères libre, mais sans blanc, pour indiquer l’endroit de la mesure
<calibration> : optionnel, à fournir si la valeur de calibration par défaut (normalement mise dans le programme) doit être remplacée par celle-ci.
“.dat” : pour indiquer un fichier de type SQM data
Remarque : les caractères “_” sont importants…

A l’entrée : tous les fichiers de mesures (type .dat), avec le nom standardisés, sont mis dans une sub-dir fixée (définition dans le programme)
A la sortie : une sub-dir fixée contiendra les “archives” des données.
Ce fichier avec nom standard : SQM_DB_<numéro de version> (en 4 positions, 0000 au départ).csv.

A chaque exécution :

  • Le dernier fichier “archive” présent est détecté et nouvelle version sera créée (0000 => 0001 => 0002, etc…)
  • le contenu de tous les fichiers détectés à l’entrée seront ajoutés au fichier “archive” courant, avec calibration et le nom de la localisation ajoutée à chaque mesure. En finale : tous les doubles éventuels sont supprimés.

Ce fichier “archive “pourra facilement être chargé dans une DB ou un programme (ex : python dataframe (Pandas)) pour des études.

# -*- coding: utf-8 -*-
"""
Created on Sat May 26 07:29:05 2018
@author: TTF

Target : add a SQM report to a global archive file for Pandas processing
"""
import pandas as pd
import glob

#initialise path & repository values, if "*" = automatic version
repositoryInputVersion  = '*'
repositoryOutputVersion = '*'
repositoryPath          = 'E:\\SQM\\$SQM_Repository\\'
inputPattern            = 'E:\\SQM\\Input\\*.dat'

#check for dynamic or static allocation for input repository
if (repositoryInputVersion == '*'): 
    repositoryInputFile=""
    for fname in glob.glob(repositoryPath+"SQM_DB_*", recursive=True):
        if fname > repositoryInputFile:
            repositoryInputFile = fname
    print("Last archive selected :",repositoryInputFile)
else:
    repositoryInputFile     = repositoryPath+'SQM_DB_'+ repositoryInputVersion +'.csv'

#check for dynamic or static allocation for output repository
if (repositoryOutputVersion == '*'): 
    elem = repositoryInputFile.split("\\")
    lastvalue=elem[len(elem)-1].split(".")[0] 
    elem = lastvalue.split("_")
    repositoryOutputVersion=str(int(elem[len(elem)-1]) + 1).zfill(4)
    print("New archive version selected :",repositoryOutputVersion)

repositoryOutputFile    = repositoryPath + "SQM_DB_" + repositoryOutputVersion +'.csv'

#open repositoryfile
sqm_DB = pd.read_csv(repositoryInputFile,delimiter=';',comment='#',
                     names=['Location','UTC','LC','Temp','Volt','Mag'])

print("Previous archive read :",sqm_DB.size,"records")

cpt_fname=0
#scan for new file
#Standard filename is <Date>"_"<time>"_"<location>("_"<calibration>)".dat"

for fname in glob.glob(inputPattern, recursive=True):
    print("Input selected :",fname)
    cpt_fname+=1
    #default calibration_value for plastic transparency (100x real value)
    calibration= 0.11
    calib_value= "0"
    
    #if present in filename, override default calibration
    try:
        path, date, time, location, calib_value = fname.split("_")
        calib_value = calib_value.split('.')[0]
    except ValueError:
        path, date, time, location = fname.split("_")
    location=location.split('.')[0]
    print("Current location : ",location)    
    
    #calibration put in filename is 100x value and must be valid
    if (float(calib_value) > 0): 
        calibration = float(calib_value) / 100 
        
    #read new SQM data file            
    df_read = pd.read_csv(fname,delimiter=';',comment='#',
                     names=['UTC','LC','Temp','Volt','Mag'])
    
    #add data to frame & add calibration on all magnitude    
    df_read['Location']=location
    df_read['Calib'] = calibration
    df_read['Mag'] = df_read['Mag'] + df_read['Calib']
    clean_df = df_read.drop(columns='Calib')
    #sort by UTC       
    df = clean_df.sort_values(by='UTC')
    #add to repository
    frames = [sqm_DB,df]
    sqm_DB = pd.concat(frames)
    print("Add ",df_read.size," records to existing",sqm_DB.size)

print("*** End add of", cpt_fname," data files added")

#suppress any duplicate (location, time, value)       
print("Before clean : ",sqm_DB.size,"records")
sqm_DB = sqm_DB.drop_duplicates(keep='first')
print("After clean  : ",sqm_DB.size,"records saved")

#save new repository
sqmout = sqm_DB.loc[:,['Location','UTC','LC','Temp','Volt','Mag']]
sqmout.to_csv(repositoryOutputFile,index=False,sep=';')
print(repositoryOutputFile,"archive saved")

A l’exécution (en considérant 3 fichiers présents à l’entrée)

Last repository selected : E:\\SQM\$SQM_Repository\SQM_DB_0000.csv
New repository value selected : 0001
Initial repository read : 0 records
Input selected : E:\SQM\Input\20180224_173752_LOC1.dat
Current location :  LOC1
Add  34181  records to existing 29298
Input selected : E:\SQM\Input\20180408_194947_LOC1.dat
Current location :  LOC1
Add  47495  records to existing 70008
Input selected : E:\SQM\Input\20200420_104312_LOC2_30.dat
Current location :  LOC2
Add  61712  records to existing 122904
*** End add of 3  data files added
Before clean :  122904 records
After clean  :  93606 records saved
E:\SQM\$SQM_Repository\SQM_DB_0001.csv repository saved

Remarque : on peut fixer le nom précis d’une version d’une archive (enlever les “*”), sinon tout est automatique par défaut…

Si on regarde le contenu, au niveau Python…
La variable “sqmOut” est un dataframe, c’est-à-dire :

Vue générale d’un ‘dataframe”

C’est une structure de table extrêmement puissante pour la librairie Pandas, qui permettra de faire des travaux de statistiques de manière simplifiée. Ainsi, une vue synthétique des données stockées s’obtient via :

sqmout.describe()
Out[3]: 
       Location                      UTC  ...   Volt    Mag
count     31203                    31203  ...  31203  31203
unique        3                    15602  ...     47   3015
top        LOC2  2018-03-19T23:00:08.000  ...   4.88    0.3
freq      17632                        2  ...   4885   4963

[4 rows x 6 columns]

Pour plus d’info, lisez ceci ou ceci