diff --git a/CHANGELOG.txt b/CHANGELOG.txt new file mode 100755 index 0000000000000000000000000000000000000000..be94ab3ff8fcc8c32c959924b1f5f922b6d98c65 --- /dev/null +++ b/CHANGELOG.txt @@ -0,0 +1,73 @@ +Version 0.4.0 + general: + port to python3 + +Version 0.3.1 + general: + Disable by default datacenter + plot: + Change default plot size + Use tight_layout + Improve detection of AM/PM dates + Allow to plot only the 2nd plot (NSB vs datetime) + plot.py now works also as a standalone tool (with user provided data file path). + + +Version 0.3.0 + general: + Added datacenter support + + +Version 0.2.2 + general: + Adopt v1.0 of the standard format + (including the filename for the daily data and plots) + read: + put the rx,cx and ix data in the header + plot: + Change de Serial number label + + +Version 0.2.1 + read: + Print the errors in make_plot call on screen. + plot: + Only print the PM/AM/Moon labels on one panel. + Print the SQM serial number. + +Version 0.2.0 + general: + Deep changes to make the program more modular. + The program now can be packaged as a single .exe file with PyInstaller. + The program can also be packaged for Linux systems. + read: + Try to use the fixed device address before looking for it automatically + this should allow the use of multiple devices in a single computer. + plot: + Code cleanup. + Use local date/time in plots. + Write statistics file. + Use pyephem to calculate the moon phase (more accurate). + Show the Moon max altitude (transit altitude or culmination). + Plot the astronomical twilights. + Object Oriented programming. + email: + Now the program can be distributed without email module. + +Version 0.1.X + read: + Variables moved to config file. + Clean-up of the code. + Improve device reset. + New read software. OO programing. + plot: + Variables moved to config file. + Renamed from plot_sqmle.py to pysqm_plot.py + Make the code and linebreaks less ugly + Fixed axis. + Moon phase plot. + email: + Renamed from email_sqmle.py to pysqm_email.py + +Version 0.0.X + First version. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000000000000000000000000000000000000..42eb4101e514e6e82ecd5111fe09e81a786de694 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1 @@ +include LICENSE.txt diff --git a/README.txt b/README.txt new file mode 100644 index 0000000000000000000000000000000000000000..ae31a9cd1b4e883910914fd54c483967e12c623b --- /dev/null +++ b/README.txt @@ -0,0 +1,115 @@ +PySQM +===== + +PySQM is a multi-platform, open-source software designed to read and +plot data from Unihedron SQM-LE and SQM-LU photometers, giving as +an output files with the 'International Dark Sky Association (IDA) +NSBM Community Standards for Reporting Skyglow Observations' format +(http://www.darksky.org/night-sky-conservation/248). + +PySQM is distributed under GNU GPL, either version 3 of the License, +or (at your option) any later version. See the file LICENSE.txt for details. + +This software has been developed by Mireia Nievas <mnievas@ucm.es> with +the invaluable help of: + + - Sergio Pascual (UCM) + - Jaime Zamorano (UCM) + - Laura Barbas (OAN) + - Pablo de Vicente (OAN) + +The initial port to Python3 has been done by Anthony Tekatch (Unihedron). + +SETUP +===== + +After downloading the software, you need to modify the file config.py. +In this file you will find several variables that need to be configured +to match your hardware settings. For example: + + - Location of the observatory (geographical coordinates). + - Device identifier. + - Device address (either IP address for SQM-LE or COM/ttyUSB port). + - Location of the data files. + - Axis limits for the plot. + +Remember that python (2.7) syntax is mandatory in this file + + +HOW TO USE THE SOFTWARE +======================= + +After configuring the software, make sure you are in the parent directory were +the README, LICENSE and MANIFEST files are located +> ls +LICENSE.txt MANIFEST.in README.txt pysqm config.py setup.py + +And then run the software. +> python -m pysqm + +The program should find your SQM device and the data adquisition.will start +(if it's night-time). + +In some systems, where python3 is the default version of python, you need +to specify python2 as the interpreter to use. This is done usually running +it as: + +> python2 -m pysqm + +or + +> python2.7 -m pysqm + +Note: running the setup.py script is neither tested nor required. +The program is currently being redesigned as a normal python package, but at +present no setup is required. + + +HOW IT WORKS +============ + +In a first step, the program tries to connect to the SQM photometer and takes +some 'tests' measures (metadata/information, calibration and data) to check +that the device is working as expected. + +After that, the program begins data acdquisition. In each iteration, it checks +whether it is night-time. In that case new data is taken. + +Each N measurements, the main program calls a plotting function to generate +a graphical representation of the current nightly data. + + +PySQM known issues +================== + +Non-ASCII characters are not supported in the config.py file. Please, avoid using 'ñ', accented vowels, etc. +In headless systems, such as the Raspberry PI, if you run the program without X, you may suffer from the following fatal error when the program tries to generate the plot: + +This application failed to start because it could not find or load the Qt platform plugin “xcb”. +Available platform plugins are: eglfs, kms, linuxfb, minimal, minimalegl, offscreen, xcb. +Reinstalling the application may fix this problem. Aborted (core dumped) + +In order to avoid this problem, you need to create (or modify if the file exists) in your HOME directory the following file: + +.config/matplotlib/matplotlibrc + +You just need to set the matplotlib backend to Agg: +backend : Agg + +Save the changes and exit. Now, PySQM should make the plots without issues. You may need to restart PySQM to apply the changes. + +Path to EXE files (windows only): +https://www.dropbox.com/s/xlbr6ktk8spjsse/PySQM.exe?dl=0 + +CHANGELOG +========= + +v0.3: + Added datacenter option (optional, disabled by default) + +v0.2: + ... + +v0.1: + ... + diff --git a/config.py b/config.py new file mode 100644 index 0000000000000000000000000000000000000000..8c4316685f5745632790ba2d4f9a520c42751b04 --- /dev/null +++ b/config.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python + +''' +PySQM configuration File. +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. + +____________________________ +Notes: + +You may need to change the following variables to match your +observatory coordinates, instrumental properties, etc. + +Python syntax is mandatory. +____________________________ +''' + + +''' +------------- +SITE location +------------- +''' + +_observatory_name = 'GURUGU' +_observatory_latitude = 40.447862 +_observatory_longitude = -3.364992 +_observatory_altitude = 680 +_observatory_horizon = 10 # If Sun is below this altitude, the program will take data + +_device_shorttype = 'SQM' # Device STR in the file +_device_type = 'SQM_LU' # Device type in the Header +_device_id = _device_type + '-' + _observatory_name # Long Device lame +_device_locationname = 'Villalbilla/Spain - Observatorio GURUGU' # Device location in the world +_data_supplier = 'Mireia Nievas / Universidad Complutense de Madrid' # Data supplier (contact) +_device_addr = '/dev/ttyUSB1' # Default IP address of the ethernet device (if not automatically found) +_measures_to_promediate = 1 # Take the mean of N measures +_delay_between_measures = 2 # Delay between two measures. In seconds. +_cache_measures = 1 # Get X measures before writing on screen/file +_plot_each = 1 # Call the plot function each X measures. + +_use_mysql = False # Set to True if you want to store data on a MySQL db. +_mysql_host = None # Host (ip:port / localhost) of the MySQL engine. +_mysql_user = None # User with write permission on the db. +_mysql_pass = None # Password for that user. +_mysql_database = None # Name of the database. +_mysql_dbtable = None # Name of the table +_mysql_port = None # Port of the MySQL server. + +_local_timezone = +1 # UTC+1 +_computer_timezone = +0 # UTC +_offset_calibration = -0.11 # magnitude = read_magnitude + offset +_reboot_on_connlost = False # Reboot if we loose connection + +# Monthly (permanent) data +monthly_data_directory = "/tmp/sqm_gurugu/" +# Daily (permanent) data +daily_data_directory = monthly_data_directory+"/datos_diarios/" +limits_nsb = [20.0,16.5] # Limits in Y-axis + +# Daily (permanent) graph +daily_graph_directory = monthly_data_directory+"/graficos_diarios/" +# Current data, deleted each day. +current_data_directory = monthly_data_directory +# Current graph, deleted each day. +current_graph_directory = monthly_data_directory +# Summary with statistics for the night +summary_data_directory = monthly_data_directory + + +''' +---------------------------- +PySQM data center (OPTIONAL) +---------------------------- +''' + +# Send the data to the data center +_send_to_datacenter = False + + +''' +Ploting options +''' +full_plot = True +limits_nsb = [20.0,16.5] # Limits in Y-axis +limits_time = [17,9] # Hours +limits_sunalt = [-80,5] # Degrees + +''' +Email options +''' +_send_data_by_email = False + diff --git a/pysqm/__init__.py b/pysqm/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..f759be3ee77037f80128d7a3062a5ba13ae3e9df --- /dev/null +++ b/pysqm/__init__.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +''' +PySQM __init__ code +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +''' + +__author__ = "Mireia Nievas" +__copyright__ = "Copyright (c) 2014 Mireia Nievas" +__credits__ = [\ + "Mireia Nievas @ UCM",\ + "Jaime Zamorano @ UCM",\ + "Laura Barbas @ OAN",\ + "Pablo de Vicente @ OAN"\ + "Anthony Tekatch @ Unihedron "\ + ] +__license__ = "GNU GPL v3" +__shortname__ = "PySQM" +__longname__ = "Python Sky Quality Meter pipeline" +__version__ = "0.4.0" +__maintainer__ = "Thorsten Alteholz" +__email__ = "python[at]alteholz[dot]de" +__status__ = "Development" # "Prototype", "Development", or "Production" + +from types import ModuleType +import sys + diff --git a/pysqm/__main__.py b/pysqm/__main__.py new file mode 100644 index 0000000000000000000000000000000000000000..1dae668549236154dbe5c9cf3e2c886988ae10b6 --- /dev/null +++ b/pysqm/__main__.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python + +''' +PySQM __main__ code +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +____________________________ +''' + +#from types import ModuleType +#import sys + +import pysqm.main as main + +while(1==1): + # Loop forever to make sure the program does not die. + try: + main.loop() + except Exception as e: + print('') + print('FATAL ERROR while running the main loop !!') + print('Error was:') + print(e) + print('Trying to restart') + print('') + diff --git a/pysqm/common.py b/pysqm/common.py new file mode 100644 index 0000000000000000000000000000000000000000..5d4a909f8ee96f53cac60294e225618e893229e6 --- /dev/null +++ b/pysqm/common.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python + +''' +PySQM common code +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +____________________________ +''' + +import math +import ephem +import datetime + +# Read the config variables from config.py +import pysqm.settings as settings +config = settings.GlobalConfig.config + +def define_ephem_observatory(): + ''' Define the Observatory in Pyephem ''' + OBS = ephem.Observer() + OBS.lat = config._observatory_latitude*ephem.pi/180 + OBS.lon = config._observatory_longitude*ephem.pi/180 + OBS.elev = config._observatory_altitude + return(OBS) + +def remove_linebreaks(data): + # Remove line breaks from data + data = data.replace('\r\n','') + data = data.replace('\r','') + data = data.replace('\n','') + return(data) + +def format_value(data,remove_str=' '): + # Remove string and spaces from data + data = remove_linebreaks(data) + data = data.replace(remove_str,'') + data = data.replace(' ','') + return(data) + +def format_value_list(data,remove_str=' '): + # Remove string and spaces from data array/list + data = [format_value(line,remove_str).split(';') for line in data] + return(data) + +def set_decimals(number,dec=3): + str_number = str(number) + int_,dec_ = str_number.split('.') + while len(dec_)<=dec: + dec_=dec_+'0' + + return(int_+'.'+dec_[:dec]) + +class observatory(object): + def read_datetime(self): + # Get UTC datetime from the computer. + utc_dt = datetime.datetime.utcnow() + #utc_dt = datetime.datetime.now() - datetime.timedelta(hours=config._computer_timezone) + #time.localtime(); daylight_saving=_.tm_isdst>0 + return(utc_dt) + + def local_datetime(self,utc_dt): + # Get Local datetime from the computer, without daylight saving. + return(utc_dt + datetime.timedelta(hours=config._local_timezone)) + + def calculate_sun_altitude(self,OBS,timeutc): + # Calculate Sun altitude + OBS.date = ephem.date(timeutc) + Sun = ephem.Sun(OBS) + return(Sun.alt) + + def next_sunset(self,OBS): + # Next sunset calculation + previous_horizon = OBS.horizon + OBS.horizon = str(config._observatory_horizon) + next_setting = OBS.next_setting(ephem.Sun()).datetime() + next_setting = next_setting.strftime("%Y-%m-%d %H:%M:%S") + OBS.horizon = previous_horizon + return(next_setting) + + def is_nighttime(self,OBS): + # Is nightime (sun below a given altitude) + timeutc = self.read_datetime() + if self.calculate_sun_altitude(OBS,timeutc)*180./math.pi>config._observatory_horizon: + return False + else: + return True + + + +RAWHeaderContent = '''# Definition of the community standard for skyglow observations 1.0 +# URL: http://www.darksky.org/NSBM/sdf1.0.pdf +# Number of header lines: 35 +# This data is released under the following license: ODbL 1.0 http://opendatacommons.org/licenses/odbl/summary/ +# Device type: $DEVICE_TYPE +# Instrument ID: $DEVICE_ID +# Data supplier: $DATA_SUPPLIER +# Location name: $LOCATION_NAME +# Position: $OBSLAT, $OBSLON, $OBSALT +# Local timezone: $TIMEZONE +# Time Synchronization: NTP +# Moving / Stationary position: STATIONARY +# Moving / Fixed look direction: FIXED +# Number of channels: 1 +# Filters per channel: HOYA CM-500 +# Measurement direction per channel: 0., 0. +# Field of view: 20 +# Number of fields per line: 6 +# SQM serial number: $SERIAL_NUMBER +# SQM firmware version: $FEATURE_NUMBER +# SQM cover offset value: $OFFSET +# SQM readout test ix: $IXREADOUT +# SQM readout test rx: $RXREADOUT +# SQM readout test cx: $CXREADOUT +# Comment: +# Comment: +# Comment: +# Comment: +# Comment: Capture program: PySQM +# blank line 30 +# blank line 31 +# blank line 32 +# UTC Date & Time, Local Date & Time, Temperature, Counts, Frequency, MSAS +# YYYY-MM-DDTHH:mm:ss.fff;YYYY-MM-DDTHH:mm:ss.fff;Celsius;number;Hz;mag/arcsec^2 +# END OF HEADER +''' diff --git a/pysqm/main.py b/pysqm/main.py new file mode 100644 index 0000000000000000000000000000000000000000..d8e2067cf3434ff00270fb77c7b733bd113c018b --- /dev/null +++ b/pysqm/main.py @@ -0,0 +1,201 @@ +#!/usr/bin/env python + +''' +PySQM main program +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +____________________________ +''' + +import os,sys +import time +import datetime +import argparse + + +''' +Read input arguments (if any) +''' +import pysqm.settings as settings +InputArguments = settings.ArgParser() +configfilename = InputArguments.config + +# Load config contents into GlobalConfig +settings.GlobalConfig.read_config_file(configfilename) + +# Get the actual config +config = settings.GlobalConfig.config + +### Load now the rest of the modules +from pysqm.read import * +import pysqm.plot + +''' +This import section is only for software build purposes. +Dont worry if some of these are missing in your setup. +''' + +def relaxed_import(themodule): + try: exec('import '+str(themodule)) + except: pass + +relaxed_import('socket') +relaxed_import('serial') +relaxed_import('_mysql') +relaxed_import('pysqm.email') + + +''' +Conditional imports +''' + +# If the old format (SQM_LE/SQM_LU) is used, replace _ with - +config._device_type = config._device_type.replace('_','-') + +if config._device_type == 'SQM-LE': + import socket +elif config._device_type == 'SQM-LU': + import serial +if config._use_mysql == True: + import _mysql + + +# Create directories if needed +for directory in [config.monthly_data_directory,config.daily_data_directory,config.current_data_directory]: + if not os.path.exists(directory): + os.makedirs(directory) + + +''' +Select the device to be used based on user input +and start the measures +''' + +if config._device_type=='SQM-LU': + mydevice = SQMLU() +elif config._device_type=='SQM-LE': + mydevice = SQMLE() +else: + print(('ERROR. Unknown device type '+str(config._device_type))) + exit(0) + + +def loop(): + ''' + Ephem is used to calculate moon position (if above horizon) + and to determine start-end times of the measures + ''' + observ = define_ephem_observatory() + niter = 0 + DaytimePrint=True + print('Starting readings ...') + while 1<2: + ''' The programs works as a daemon ''' + utcdt = mydevice.read_datetime() + #print (str(mydevice.local_datetime(utcdt))), + if mydevice.is_nighttime(observ): + # If we are in a new night, create the new file. + config._send_to_datacenter = False ### Not enabled by default + try: + assert(config._send_to_datacenter == True) + assert(niter == 0) + mydevice.save_data_datacenter("NEWFILE") + except: pass + + StartDateTime = datetime.datetime.now() + niter += 1 + + mydevice.define_filenames() + + ''' Get values from the photometer ''' + try: + timeutc_mean,timelocal_mean,temp_sensor,\ + freq_sensor,ticks_uC,sky_brightness = \ + mydevice.read_photometer(\ + Nmeasures=config._measures_to_promediate,PauseMeasures=10) + except: + print('Connection lost') + if config._reboot_on_connlost == True: + sleep(600) + os.system('reboot.bat') + + time.sleep(1) + mydevice.reset_device() + + formatted_data = mydevice.format_content(\ + timeutc_mean,timelocal_mean,temp_sensor,\ + freq_sensor,ticks_uC,sky_brightness) + + try: + assert(config._use_mysql == True) + mydevice.save_data_mysql(formatted_data) + except: pass + + try: + assert(config._send_to_datacenter == True) + mydevice.save_data_datacenter(formatted_data) + except: pass + + mydevice.data_cache(formatted_data,number_measures=config._cache_measures,niter=niter) + + if niter%config._plot_each == 0: + ''' Each X minutes, plot a new graph ''' + try: pysqm.plot.make_plot(send_emails=False,write_stats=False) + except: + print('Warning: Error plotting data.') + print((sys.exc_info())) + + if DaytimePrint==False: + DaytimePrint=True + + MainDeltaSeconds = (datetime.datetime.now()-StartDateTime).total_seconds() + time.sleep(max(1,config._delay_between_measures-MainDeltaSeconds)) + + else: + ''' Daytime, print info ''' + if DaytimePrint==True: + utcdt = utcdt.strftime("%Y-%m-%d %H:%M:%S") + print((utcdt), end=' ') + print(('. Daytime. Waiting until '+str(mydevice.next_sunset(observ)))) + DaytimePrint=False + if niter>0: + mydevice.flush_cache() + if config._send_data_by_email==True: + try: pysqm.plot.make_plot(send_emails=True,write_stats=True) + except: + print('Warning: Error plotting data / sending email.') + print((sys.exc_info())) + + else: + try: pysqm.plot.make_plot(send_emails=False,write_stats=True) + except: + print('Warning: Error plotting data.') + print((sys.exc_info())) + + niter = 0 + + # Send data that is still in the datacenter buffer + try: + assert(config._send_to_datacenter == True) + mydevice.save_data_datacenter("") + except: pass + + time.sleep(300) + + diff --git a/pysqm/plot.py b/pysqm/plot.py new file mode 100644 index 0000000000000000000000000000000000000000..8f127a23ec8f4aa7604f8d0b965ac7a232a28d9a --- /dev/null +++ b/pysqm/plot.py @@ -0,0 +1,879 @@ +#!/usr/bin/env python + +''' +PySQM plotting program +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +____________________________ +''' + +import os,sys +import ephem +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.ticker as ticker +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +from datetime import datetime,timedelta + +# Read configuration + +if __name__ != '__main__': + import pysqm.settings as settings + config = settings.GlobalConfig.config + + for directory in [\ + config.monthly_data_directory,\ + config.daily_graph_directory,\ + config.current_graph_directory]: + if not os.path.exists(directory): + os.makedirs(directory) + + +class Ephemerids(object): + def __init__(self): + from pysqm.common import define_ephem_observatory + self.Observatory = define_ephem_observatory() + + def ephem_date_to_datetime(self,ephem_date): + # Convert ephem dates to datetime + date_,time_ = str(ephem_date).split(' ') + date_ = date_.split('/') + time_ = time_.split(':') + + return(datetime(\ + int(date_[0]),int(date_[1]),int(date_[2]),\ + int(time_[0]),int(time_[1]),int(time_[2]))) + + def end_of_the_day(self,thedate): + newdate = thedate+timedelta(days=1) + newdatetime = datetime(\ + newdate.year,\ + newdate.month,\ + newdate.day,0,0,0) + newdatetime = newdatetime-timedelta(hours=config._local_timezone) + return(newdatetime) + + + def calculate_moon_ephems(self,thedate): + # Moon ephemerids + self.Observatory.horizon = '0' + self.Observatory.date = str(self.end_of_the_day(thedate)) + + # Moon phase + Moon = ephem.Moon() + Moon.compute(self.Observatory) + self.moon_phase = Moon.phase + self.moon_maxelev = Moon.transit_alt + + try: + float(self.moon_maxelev) + except: + # The moon has no culmination time for 1 day + # per month, so there is no max altitude. + # As a workaround, we use the previous day culmination. + # The error should be small. + + # Set the previous day date + thedate2 = thedate - timedelta(days=1) + self.Observatory.date = str(self.end_of_the_day(thedate2)) + Moon2 = ephem.Moon() + Moon2.compute(self.Observatory) + self.moon_maxelev = Moon2.transit_alt + + # Recover the real date + self.Observatory.date = str(self.end_of_the_day(thedate)) + + # Moon rise and set + self.moon_prev_rise = \ + self.ephem_date_to_datetime(self.Observatory.previous_rising(ephem.Moon())) + self.moon_prev_set = \ + self.ephem_date_to_datetime(self.Observatory.previous_setting(ephem.Moon())) + self.moon_next_rise = \ + self.ephem_date_to_datetime(self.Observatory.next_rising(ephem.Moon())) + self.moon_next_set = \ + self.ephem_date_to_datetime(self.Observatory.next_setting(ephem.Moon())) + + def calculate_twilight(self,thedate,twilight=-18): + ''' + Changing the horizon forces ephem to + calculate different types of twilights: + -6: civil, + -12: nautical, + -18: astronomical, + ''' + self.Observatory.horizon = str(twilight) + self.Observatory.date = str(self.end_of_the_day(thedate)) + + self.twilight_prev_rise = self.ephem_date_to_datetime(\ + self.Observatory.previous_rising(ephem.Sun(),use_center=True)) + self.twilight_prev_set = self.ephem_date_to_datetime(\ + self.Observatory.previous_setting(ephem.Sun(),use_center=True)) + self.twilight_next_rise = self.ephem_date_to_datetime(\ + self.Observatory.next_rising(ephem.Sun(),use_center=True)) + self.twilight_next_set = self.ephem_date_to_datetime(\ + self.Observatory.next_setting(ephem.Sun(),use_center=True)) + + + +class SQMData(object): + # Split pre and after-midnight data + + class premidnight(object): + pass + + class aftermidnight(object): + pass + + class Statistics(object): + pass + + def __init__(self,filename,Ephem): + self.all_night_sb = [] + self.all_night_dt = [] + self.all_night_temp = [] + + for variable in [\ + 'utcdates','localdates','sun_altitudes',\ + 'temperatures','tick_counts','frequencies',\ + 'night_sbs','label_dates','sun_altitude']: + setattr(self.premidnight,variable,[]) + setattr(self.aftermidnight,variable,[]) + + self.load_rawdata(filename) + self.process_rawdata(Ephem) + self.check_number_of_nights() + + def extract_metadata(self,raw_data_and_metadata): + from pysqm.common import format_value + metadata_lines = [\ + line for line in raw_data_and_metadata \ + if format_value(line)[0]=='#'] + + # Extract the serial number + serial_number_line = [\ + line for line in metadata_lines \ + if 'SQM serial number:' in line][0] + self.serial_number = format_value(serial_number_line.split(':')[-1]) + + def check_validdata(self,data_line): + from pysqm.common import format_value + try: + assert(format_value(data_line)[0]!='#') + assert(format_value(data_line)[0]!='') + except: + return(False) + else: + return(True) + + def load_rawdata(self,filename): + ''' + Open the file, read the data and close the file + ''' + sqm_file = open(filename, 'r') + raw_data_and_metadata = sqm_file.readlines() + self.metadata = self.extract_metadata(raw_data_and_metadata) + + self.raw_data = [\ + line for line in raw_data_and_metadata \ + if self.check_validdata(line)==True] + sqm_file.close() + + def process_datetimes(self,str_datetime): + ''' + Get date and time in a str format + Return as datetime object + ''' + str_date,str_time = str_datetime.split('T') + + year = int(str_date.split('-')[0]) + month = int(str_date.split('-')[1]) + day = int(str_date.split('-')[2]) + + # Time may be not complete. Workaround + hour = int(str_time.split(':')[0]) + try: + minute = int(str_time.split(':')[1]) + except: + minute = 0 + second = 0 + else: + try: + second = int(str_time.split(':')[2]) + except: + second = 0 + + return(datetime(year,month,day,hour,minute,second)) + + def process_rawdata(self,Ephem): + from pysqm.common import format_value_list + ''' + Get the important information from the raw_data + and put it in a more useful format + ''' + self.raw_data = format_value_list(self.raw_data) + + for k,line in enumerate(self.raw_data): + # DateTime extraction + utcdatetime = self.process_datetimes(line[0]) + localdatetime = self.process_datetimes(line[1]) + + # Check that datetimes are corrent + calc_localdatetime = utcdatetime+timedelta(hours=config._local_timezone) + if (calc_localdatetime != localdatetime): return 1 + + # Set the datetime for astronomical calculations. + Ephem.Observatory.date = ephem.date(utcdatetime) + + # Date in str format: 20130115 + label_date = str(localdatetime.date()).replace('-','') + + # Temperature + temperature = float(line[2]) + # Counts + tick_counts = float(line[3]) + # Frequency + frequency = float(line[4]) + # Night sky background + night_sb = float(line[5]) + try: config._plot_corrected_nsb + except AttributeError: config._plot_corrected_data=False + if (config._plot_corrected_data): + night_sb += config._plot_corrected_data*config._offset_calibration + # Define sun in pyephem + Sun = ephem.Sun(Ephem.Observatory) + + self.premidnight.label_date=[] + self.aftermidnight.label_dates=[] + + + if localdatetime.hour > 12: + self.premidnight.utcdates.append(utcdatetime) + self.premidnight.localdates.append(localdatetime) + self.premidnight.temperatures.append(temperature) + self.premidnight.tick_counts.append(tick_counts) + self.premidnight.frequencies.append(frequency) + self.premidnight.night_sbs.append(night_sb) + self.premidnight.sun_altitude.append(Sun.alt) + if label_date not in self.premidnight.label_dates: + self.premidnight.label_dates.append(label_date) + else: + self.aftermidnight.utcdates.append(utcdatetime) + self.aftermidnight.localdates.append(localdatetime) + self.aftermidnight.temperatures.append(temperature) + self.aftermidnight.tick_counts.append(tick_counts) + self.aftermidnight.frequencies.append(frequency) + self.aftermidnight.night_sbs.append(night_sb) + self.aftermidnight.sun_altitude.append(Sun.alt) + if label_date not in self.aftermidnight.label_dates: + self.aftermidnight.label_dates.append(label_date) + + # Data for the complete night + self.all_night_dt.append(utcdatetime) # Must be in UTC! + self.all_night_sb.append(night_sb) + self.all_night_temp.append(temperature) + + def check_number_of_nights(self): + ''' + Check that the number of nights is exactly 1 and + extract it to a new variable self.Night. + Needed for the statistics part of the analysis and + to make the plot. + ''' + + if np.size(self.premidnight.localdates)>0: + self.Night = np.unique([DT.date() \ + for DT in self.premidnight.localdates])[0] + elif np.size(self.aftermidnight.localdates)>0: + self.Night = np.unique([(DT-timedelta(hours=12)).date() \ + for DT in self.aftermidnight.localdates])[0] + else: + print('Warning, No Night detected.') + self.Night = None + + def data_statistics(self,Ephem): + ''' + Make statistics on the data. + Useful to summarize night conditions. + ''' + def select_bests(values,number): + return(np.sort(values)[::-1][0:number]) + + def fourier_filter(array,nterms): + ''' + Make a fourier filter for the first nterms terms. + ''' + array_fft = np.fft.fft(array) + # Filter data + array_fft[nterms:]=0 + filtered_array = np.fft.ifft(array_fft) + return(filtered_array) + + def window_smooth(x,window_len=10,window='hanning'): + # http://scipy-cookbook.readthedocs.io/items/SignalSmooth.html + if x.ndim != 1: raise ValueError("smooth requires 1-d arrays") + if x.size < window_len: raise ValueError("size(input) < window_size") + if window_len < 3: return x + if not window in ['flat','hanning','hamming','bartlett','blackman']: + raise ValueError("Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'") + s=np.r_[x[window_len-1:0:-1],x,x[-2:-window_len-1:-1]] + if window == 'flat': #moving average + w=np.ones(window_len,'d') + else: + w=eval('np.'+window+'(window_len)') + y=np.convolve(w/w.sum(),s,mode='valid') + return(y) + + astronomical_night_filter = (\ + (np.array(self.all_night_dt)>Ephem.twilight_prev_set)*\ + (np.array(self.all_night_dt)<Ephem.twilight_next_rise)) + + if np.sum(astronomical_night_filter)>10: + self.astronomical_night_sb = \ + np.array(self.all_night_sb)[astronomical_night_filter] + self.astronomical_night_temp = \ + np.array(self.all_night_temp)[astronomical_night_filter] + else: + print((\ + 'Warning, < 10 points in astronomical night, '+\ + ' using the whole night data instead')) + self.astronomical_night_sb = self.all_night_sb + self.astronomical_night_temp = self.all_night_temp + + Stat = self.Statistics + #with self.Statistics as Stat: + # Complete list + Stat.mean = np.mean(self.astronomical_night_sb) + Stat.median = np.median(self.astronomical_night_sb) + Stat.std = np.median(self.astronomical_night_sb) + Stat.number = np.size(self.astronomical_night_sb) + # Only the best 1/100th. + Stat.bests_number = round(1+Stat.number/25) + Stat.bests_mean = np.mean(select_bests(self.astronomical_night_sb,Stat.bests_number)) + Stat.bests_median = np.median(select_bests(self.astronomical_night_sb,Stat.bests_number)) + Stat.bests_std = np.std(select_bests(self.astronomical_night_sb,Stat.bests_number)) + Stat.bests_err = Stat.bests_std*1./np.sqrt(Stat.bests_number) + + Stat.model_nterm = round(1+Stat.number/25) + #data_smooth = fourier_filter(self.astronomical_night_sb,nterms=Stat.model_nterm) + data_smooth = window_smooth(self.astronomical_night_sb, + window_len=Stat.model_nterm) + min_length = min(len(data_smooth),len(self.astronomical_night_sb)) + data_residuals = self.astronomical_night_sb[:min_length]-data_smooth[:min_length] + Stat.data_model_abs_meandiff = np.mean(np.abs(data_residuals)) + Stat.data_model_sum_squareresiduals = np.sum(data_residuals**2) + + # Other interesting data + Stat.min_temperature = np.min(self.astronomical_night_temp) + Stat.max_temperature = np.max(self.astronomical_night_temp) + + +class Plot(object): + def __init__(self,Data,Ephem): + plt.clf() # plt.hold(True/False) is deprecated + Data = self.prepare_plot(Data,Ephem) + + try: config.full_plot + except: config.full_plot = False + if (config.full_plot): + self.make_figure(thegraph_altsun=True,thegraph_time=True) + self.plot_data_sunalt(Data,Ephem) + else: + self.make_figure(thegraph_altsun=False,thegraph_time=True) + + self.plot_data_time(Data,Ephem) + + self.plot_moonphase(Ephem) + self.plot_twilight(Ephem) + + def plot_moonphase(self,Ephem): + ''' + shade the period of time for which the moon is above the horizon + ''' + if Ephem.moon_next_rise > Ephem.moon_next_set: + # We need to divide the plotting in two phases + #(pre-midnight and after-midnight) + self.thegraph_time.axvspan(\ + Ephem.moon_prev_rise+timedelta(hours=config._local_timezone),\ + Ephem.moon_next_set+timedelta(hours=config._local_timezone),\ + edgecolor='#d62728',facecolor='#d62728', alpha=0.1,clip_on=True) + else: + self.thegraph_time.axvspan(\ + Ephem.moon_prev_rise+timedelta(hours=config._local_timezone),\ + Ephem.moon_prev_set+timedelta(hours=config._local_timezone),\ + edgecolor='#d62728',facecolor='#d62728', alpha=0.1,clip_on=True) + self.thegraph_time.axvspan(\ + Ephem.moon_next_rise+timedelta(hours=config._local_timezone),\ + Ephem.moon_next_set+timedelta(hours=config._local_timezone),\ + edgecolor='#d62728',facecolor='#d62728', alpha=0.1,clip_on=True) + + def plot_twilight(self,Ephem): + ''' + Plot vertical lines on the astronomical twilights + ''' + self.thegraph_time.axvline(\ + Ephem.twilight_prev_set+timedelta(hours=config._local_timezone),\ + color='black', ls='dashdot', lw=1, alpha=0.75, clip_on=True) + self.thegraph_time.axvline(\ + Ephem.twilight_next_rise+timedelta(hours=config._local_timezone),\ + color='black', ls='dashdot', lw=1, alpha=0.75, clip_on=True) + + def make_subplot_sunalt(self,twinplot=0): + ''' + Make a subplot. + If twinplot = 0, then this will be the only plot in the figure + if twinplot = 1, this will be the first subplot + if twinplot = 2, this will be the second subplot + ''' + if twinplot == 0: + self.thegraph_sunalt = self.thefigure.add_subplot(1,1,1) + else: + self.thegraph_sunalt = self.thefigure.add_subplot(2,1,twinplot) + + + self.thegraph_sunalt.set_title(\ + 'Sky Brightness ('+config._device_shorttype+'-'+\ + config._observatory_name+')\n',fontsize='x-large') + self.thegraph_sunalt.set_xlabel('Solar altitude (deg)',fontsize='large') + self.thegraph_sunalt.set_ylabel('Sky Brightness (mag/arcsec2)',fontsize='medium') + + # Auxiliary plot (Temperature) + ''' + self.thegraph_sunalt_temp = self.thegraph_sunalt.twinx() + self.thegraph_sunalt_temp.set_ylim(-10, 50) + self.thegraph_sunalt_temp.set_ylabel('Temperature (C)',fontsize='medium') + ''' + + # format the ticks (frente a alt sol) + tick_values = list(range(config.limits_sunalt[0],config.limits_sunalt[1]+5,5)) + tick_marks = np.multiply([deg for deg in tick_values],np.pi/180.0) + tick_labels = [str(deg) for deg in tick_values] + + self.thegraph_sunalt.set_xticks(tick_marks) + self.thegraph_sunalt.set_xticklabels(tick_labels) + self.thegraph_sunalt.yaxis.set_minor_locator(ticker.MultipleLocator(0.5)) + self.thegraph_sunalt.grid(True,which='major', + alpha=0.2,color='k',ls='',lw=0.5) + self.thegraph_sunalt.grid(True,which='minor', + alpha=0.2,color='k',ls='solid',lw=0.5) + + def make_subplot_time(self,twinplot=0): + ''' + Make a subplot. + If twinplot = 0, then this will be the only plot in the figure + if twinplot = 1, this will be the first subplot + if twinplot = 2, this will be the second subplot + ''' + if twinplot == 0: + self.thegraph_time = self.thefigure.add_subplot(1,1,1) + else: + self.thegraph_time = self.thefigure.add_subplot(2,1,twinplot) + + if config._local_timezone<0: + UTC_offset_label = '-'+str(abs(config._local_timezone)) + elif config._local_timezone>0: + UTC_offset_label = '+'+str(abs(config._local_timezone)) + else: UTC_offset_label = '' + + #self.thegraph_time.set_title('Sky Brightness (SQM-'+config._observatory_name+')',\ + # fontsize='x-large') + self.thegraph_time.set_xlabel('Time (UTC'+UTC_offset_label+')',fontsize='large') + self.thegraph_time.set_ylabel('Sky Brightness (mag/arcsec2)',fontsize='medium') + + # Auxiliary plot (Temperature) + ''' + self.thegraph_time_temp = self.thegraph_time.twinx() + self.thegraph_time_temp.set_ylim(-10, 50) + self.thegraph_time_temp.set_ylabel('Temperature (C)',fontsize='medium') + ''' + + # format the ticks (vs time) + daylocator = mdates.HourLocator(byhour=[4,20]) + hourlocator = mdates.HourLocator() + dayFmt = mdates.DateFormatter('%d %b %Y') + hourFmt = mdates.DateFormatter('%H') + + self.thegraph_time.xaxis.set_major_locator(daylocator) + self.thegraph_time.xaxis.set_major_formatter(dayFmt) + self.thegraph_time.xaxis.set_minor_locator(hourlocator) + self.thegraph_time.xaxis.set_minor_formatter(hourFmt) + self.thegraph_time.yaxis.set_minor_locator(ticker.MultipleLocator(0.5)) + + self.thegraph_time.xaxis.set_tick_params(which='major', pad=15) + self.thegraph_time.format_xdata = mdates.DateFormatter('%Y-%m-%d_%H:%M:%S') + self.thegraph_time.grid(True,which='major', + alpha=0.2,color='k',ls='',lw=0.5) + self.thegraph_time.grid(True,which='minor', + alpha=0.2,color='k',ls='solid',lw=0.5) + + def make_figure(self,thegraph_altsun=True,thegraph_time=True): + # Make the figure and the graph + if thegraph_time==False: + self.thefigure = plt.figure(figsize=(7,3.)) + self.make_subplot_sunalt(twinplot=0) + elif thegraph_altsun==False: + self.thefigure = plt.figure(figsize=(7,3.)) + self.make_subplot_time(twinplot=0) + else: + self.thefigure = plt.figure(figsize=(7,6.)) + self.make_subplot_sunalt(twinplot=1) + self.make_subplot_time(twinplot=2) + + # Adjust the space between plots + plt.subplots_adjust(hspace=0.35) + + + def prepare_plot(self,Data,Ephem): + ''' + Warning! Multiple night plot implementation is pending. + Until the support is implemented, check that no more than 1 night + is used + ''' + + # Mean datetime + dts = Data.all_night_dt + mean_dt = dts[0]+np.sum(np.array(dts)-dts[0])/np.size(dts) + sel_night = (mean_dt - timedelta(hours=12)).date() + + Data.premidnight.filter = np.array(\ + [Date.date()==sel_night for Date in Data.premidnight.localdates]) + Data.aftermidnight.filter = np.array(\ + [(Date-timedelta(days=1)).date()==sel_night\ + for Date in Data.aftermidnight.localdates]) + + return(Data) + + def plot_data_sunalt(self,Data,Ephem): + ''' + Plot NSB data vs Sun altitude + ''' + # Plot the data + TheData = Data.premidnight + if np.size(TheData.filter)>0: + self.thegraph_sunalt.plot(\ + np.array(TheData.sun_altitude)[TheData.filter],\ + np.array(TheData.night_sbs)[TheData.filter],color='#2ca02c') + ''' + self.thegraph_sunalt.plot(\ + np.array(TheData.sun_altitude)[TheData.filter],\ + np.array(TheData.temperatures)[TheData.filter],color='#9467bd',alpha=0.5)) + ''' + TheData = Data.aftermidnight + if np.size(TheData.filter)>0: + self.thegraph_sunalt.plot(\ + np.array(TheData.sun_altitude)[TheData.filter],\ + np.array(TheData.night_sbs)[TheData.filter],color='#1f77b4') + ''' + self.thegraph_sunalt.plot(\ + np.array(TheData.sun_altitude)[TheData.filter],\ + np.array(TheData.temperatures)[TheData.filter],color='#9467bd',alpha=0.5)) + ''' + + # Make limits on data range. + self.thegraph_sunalt.set_xlim([\ + config.limits_sunalt[0]*np.pi/180.,\ + config.limits_sunalt[1]*np.pi/180.]) + self.thegraph_sunalt.set_ylim(config.limits_nsb) + + premidnight_label = str(Data.premidnight.label_dates).replace('[','').replace(']','') + aftermidnight_label = str(Data.aftermidnight.label_dates).replace('[','').replace(']','') + + self.thegraph_sunalt.text(0.00,1.015,\ + config._device_shorttype+'-'+config._observatory_name+' '*5+'Serial #'+str(Data.serial_number),\ + color='0.25',fontsize='small',fontname='monospace',\ + transform = self.thegraph_sunalt.transAxes) + + self.thegraph_sunalt.text(0.75,0.92,'PM: '+premidnight_label,\ + color='#2ca02c',fontsize='small',transform = self.thegraph_sunalt.transAxes) + self.thegraph_sunalt.text(0.75,0.84,'AM: '+aftermidnight_label,\ + color='#1f77b4',fontsize='small',transform = self.thegraph_sunalt.transAxes) + + ''' + if np.size(Data.Night)==1: + self.thegraph_sunalt.text(0.75,1.015,'Moon: %d%s (%d%s)' \ + %(Ephem.moon_phase, "%", Ephem.moon_maxelev*180./np.pi,"$^\mathbf{o}$"),\ + color='#d62728',fontsize='small',fontname='monospace',\ + transform = self.thegraph_sunalt.transAxes) + ''' + + def plot_data_time(self,Data,Ephem): + ''' + Plot NSB data vs Sun altitude + ''' + + # Plot the data (NSB and temperature) + TheData = Data.premidnight + if np.size(TheData.filter)>0: + self.thegraph_time.plot(\ + np.array(TheData.localdates)[TheData.filter],\ + np.array(TheData.night_sbs)[TheData.filter],color='#2ca02c') + ''' + self.thegraph_time_temp.plot(\ + np.array(TheData.localdates)[TheData.filter],\ + np.array(TheData.temperatures)[TheData.filter],color='#9467bd',alpha=0.5) + ''' + + + TheData = Data.aftermidnight + if np.size(TheData.filter)>0: + self.thegraph_time.plot(\ + np.array(TheData.localdates)[TheData.filter],\ + np.array(TheData.night_sbs)[TheData.filter],color='#1f77b4') + ''' + self.thegraph_time_temp.plot(\ + np.array(TheData.localdates)[TheData.filter],\ + np.array(TheData.temperatures)[TheData.filter],color='#9467bd',alpha=0.5) + ''' + + # Vertical line to mark 0h + self.thegraph_time.axvline(\ + Data.Night+timedelta(days=1), + color='black', alpha=0.75,lw=1,ls='solid',clip_on=True) + + # Set the xlimit for the time plot. + if np.size(Data.premidnight.filter)>0: + begin_plot_dt = Data.premidnight.localdates[-1] + begin_plot_dt = datetime(\ + begin_plot_dt.year,\ + begin_plot_dt.month,\ + begin_plot_dt.day,\ + config.limits_time[0],0,0) + end_plot_dt = begin_plot_dt+timedelta(\ + hours=24+config.limits_time[1]-config.limits_time[0]) + elif np.size(Data.aftermidnight.filter)>0: + end_plot_dt = Data.aftermidnight.localdates[-1] + end_plot_dt = datetime(\ + end_plot_dt.year,\ + end_plot_dt.month,\ + end_plot_dt.day,\ + config.limits_time[1],0,0) + begin_plot_dt = end_plot_dt-timedelta(\ + hours=24+config.limits_time[1]-config.limits_time[0]) + else: + print('Warning: Cannot calculate plot limits') + return(None) + + self.thegraph_time.set_xlim(begin_plot_dt,end_plot_dt) + self.thegraph_time.set_ylim(config.limits_nsb) + + premidnight_label = str(Data.premidnight.label_dates).replace('[','').replace(']','') + aftermidnight_label = str(Data.aftermidnight.label_dates).replace('[','').replace(']','') + + self.thegraph_time.text(0.00,1.015,\ + config._device_shorttype+'-'+config._observatory_name+' '*5+'Serial #'+str(Data.serial_number),\ + color='0.25',fontsize='small',fontname='monospace',\ + transform = self.thegraph_time.transAxes) + + if np.size(Data.Night)==1: + self.thegraph_time.text(0.75,1.015,'Moon: %d%s (%d%s)' \ + %(Ephem.moon_phase, "%", Ephem.moon_maxelev*180./np.pi,"$^\mathbf{o}$"),\ + color='black',fontsize='small',fontname='monospace',\ + transform = self.thegraph_time.transAxes) + + def save_figure(self,output_filename): + self.thefigure.savefig(output_filename, bbox_inches='tight',dpi=150) + + def show_figure(self): + plt.show(self.thefigure) + + def close_figure(self): + plt.close('all') + + +def save_stats_to_file(Night,NSBData,Ephem): + from pysqm.common import set_decimals + ''' + Save statistics to file + ''' + + Stat = NSBData.Statistics + + Header = \ + '# Summary statistics for '+str(config._device_shorttype+'_'+config._observatory_name)+'\n'+\ + '# Description of columns (CSV file):\n'+\ + '# Col 1: Date\n'+\ + '# Col 2: Total measures\n'+\ + '# Col 3: Number of Best NSB measures\n'+\ + '# Col 4: Median of best N NSBs (mag/arcsec2)\n'+\ + '# Col 5: Err in the median of best N NSBs (mag/arcsec2)\n'+\ + '# Col 6: Window size for the smoothing function\n'+\ + '# Col 7: Mean of Abs diff of NSBs data - fourier model (mag/arcsec2)\n'+\ + '# Col 8: Min Temp (C) between astronomical twilights\n'+\ + '# Col 9: Max Temp (C) between astronomical twilights\n\n' + #'# Col 6: Number of terms of the low-freq fourier model\n'+\ + + formatted_data = \ + str(Night)+';'+\ + str(Stat.number)+';'+\ + str(Stat.bests_number)+';'+\ + set_decimals(Stat.bests_median,4)+';'+\ + set_decimals(Stat.bests_err,4)+';'+\ + str(Stat.model_nterm)+';'+\ + set_decimals(Stat.data_model_abs_meandiff,4)+';'+\ + set_decimals(Stat.min_temperature,1)+';'+\ + set_decimals(Stat.max_temperature,1)+\ + '\n' + + statistics_filename = \ + config.summary_data_directory+'/Statistics_'+\ + str(config._device_shorttype+'_'+config._observatory_name)+'.dat' + + print('Writing statistics file') + + def safe_create_file(filename): + if not os.path.exists(filename): + open(filename, 'w').close() + + def read_file(filename): + thefile = open(filename,'r') + content = thefile.read() + thefile.close() + return(content) + + def write_file(filename,content): + thefile = open(filename,'w') + thefile.write(content) + thefile.close() + + def append_file(filename,content): + thefile = open(filename,'a') + thefile.write(content) + thefile.close() + + # Create file if not exists + safe_create_file(statistics_filename) + + # Read the content + stat_file_content = read_file(statistics_filename) + + # If the file doesnt have a proper header, add it to the beginning + + def valid_line(line): + if '#' in line: + return False + elif line.replace(' ','')=='': + return False + else: + return True + + if Header not in stat_file_content: + stat_file_content = [line for line in stat_file_content.split('\n') \ + if valid_line(line)] + stat_file_content = '\n'.join(stat_file_content) + stat_file_content = Header+stat_file_content + write_file(statistics_filename,stat_file_content) + + # Remove any previous statistic for the given Night in the file + if str(Night) in stat_file_content: + stat_file_content = [line for line in stat_file_content.split('\n') \ + if str(Night) not in line] + stat_file_content = '\n'.join(stat_file_content) + write_file(statistics_filename,stat_file_content) + + # Append to the end of the file + append_file(statistics_filename,formatted_data) + + +def make_plot(input_filename=None,send_emails=False,write_stats=False): + ''' + Main function (allows to execute the program + from within python. + - Extracts the NSB data from a given data file + - Performs statistics + - Save statistics to file + - Create the plot + ''' + + print('Plotting photometer data ...') + + if (input_filename is None): + input_filename = config.current_data_directory+\ + '/'+config._device_shorttype+'_'+config._observatory_name+'.dat' + + # Define the observatory in ephem + Ephem = Ephemerids() + + # Get and process the data from input_filename + NSBData = SQMData(input_filename,Ephem) + + # Moon and twilight ephemerids. + Ephem.calculate_moon_ephems(thedate=NSBData.Night) + Ephem.calculate_twilight(thedate=NSBData.Night) + + # Calculate data statistics + NSBData.data_statistics(Ephem) + + # Write statiscs to file? + if write_stats==True: + save_stats_to_file(NSBData.Night,NSBData,Ephem) + + # Plot the data and save the resulting figure + NSBPlot = Plot(NSBData,Ephem) + + output_filenames = [\ + str("%s/%s_%s.png" %(config.current_data_directory,\ + config._device_shorttype,config._observatory_name)),\ + str("%s/%s_120000_%s-%s.png" \ + %(config.daily_graph_directory, str(NSBData.Night).replace('-',''),\ + config._device_shorttype, config._observatory_name))\ + ] + + for output_filename in output_filenames: + NSBPlot.save_figure(output_filename) + + # Close figure + NSBPlot.close_figure() + + if send_emails == True: + import pysqm.email + night_label = str(datetime.date.today()-timedelta(days=1)) + pysqm.email.send_emails(night_label=night_label,Stat=NSBData.Statistics) + +''' +The following code allows to execute plot.py as a standalone program. +''' +if __name__ == '__main__': + # Exec the main program + import settings as settings + InputArguments = settings.ArgParser(inputfile=True) + configfilename = InputArguments.config + try: + settings.GlobalConfig.read_config_file(configfilename) + config = settings.GlobalConfig.config + make_plot(input_filename=InputArguments.input,\ + send_emails=False,write_stats=True) + except: + raise + print("Error: The arguments you provided are invalid") + InputArguments.print_help() + + + + + + + + + + + + + + diff --git a/pysqm/read.py b/pysqm/read.py new file mode 100644 index 0000000000000000000000000000000000000000..ac26a1425d7113672fe8792fa71aae46176c4721 --- /dev/null +++ b/pysqm/read.py @@ -0,0 +1,854 @@ + +#!/usr/bin/env python + +''' +PySQM reading program +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +____________________________ +''' + +import os,sys +import inspect +import time +import datetime +import numpy as np +import struct +import socket + +# Default, to ignore the length of the read string. +_cal_len_ = None +_meta_len_ = None +_data_len_ = None + +from pysqm.common import * + +''' +This import section is only for software build purposes. +Dont worry if some of these are missing in your setup. +''' + +def relaxed_import(themodule): + try: exec('import '+str(themodule)) + except: pass + +relaxed_import('serial') +relaxed_import('_mysql') +relaxed_import('pysqm.email') + +''' +Read configuration +''' +import pysqm.settings as settings +config = settings.GlobalConfig.config + +try: + DEBUG=config.DEBUG +except: + DEBUG=False + +''' +Conditional imports +''' + +# If the old format (SQM_LE/SQM_LU) is used, replace _ with - +config._device_type = config._device_type.replace('_','-') + +if config._device_type == 'SQM-LE': + import socket +elif config._device_type == 'SQM-LU': + import serial +if config._use_mysql == True: + import _mysql + + +def filtered_mean(array,sigma=3): + # Our data probably contains outliers, filter them + # Notes: + # Median is more robust than mean + # Std increases if the outliers are far away from real values. + # We need to limit the amount of discrepancy we want in the data (20%?). + + # We will use data masking and some operations with arrays. Convert to numpy. + array = np.array(array) + + # Get the median and std. + data_median = np.median(array) + data_std = np.std(array) + + # Max discrepancy we allow. + fixed_max_dev = 0.2*data_median + clip_deviation = np.min([fixed_max_dev,data_std*sigma+0.1]) + + # Create the filter (10% flux + variable factor) + filter_values_ok = np.abs(array-data_median)<=clip_deviation + filtered_values = array[filter_values_ok] + + # Return the mean of filtered data or the median. + if np.size(filtered_values)==0: + print('Warning: High dispersion found on last measures') + filtered_mean = data_median + else: + filtered_mean = np.mean(filtered_values) + + return(filtered_mean) + + +class device(observatory): + def standard_file_header(self): + # Data Header, at the end of this script. + header_content=RAWHeaderContent + + # Update data file header with observatory data + header_content = header_content.replace(\ + '$DEVICE_TYPE',str(config._device_type)) + header_content = header_content.replace(\ + '$DEVICE_ID',str(config._device_id)) + header_content = header_content.replace(\ + '$DATA_SUPPLIER',str(config._data_supplier)) + header_content = header_content.replace(\ + '$LOCATION_NAME',str(config._device_locationname)) + header_content = header_content.replace(\ + '$OBSLAT',str(config._observatory_latitude)) + header_content = header_content.replace(\ + '$OBSLON',str(config._observatory_longitude)) + header_content = header_content.replace(\ + '$OBSALT',str(config._observatory_altitude)) + header_content = header_content.replace(\ + '$OFFSET',str(config._offset_calibration)) + + if config._local_timezone==0: + header_content = header_content.replace(\ + '$TIMEZONE','UTC') + elif config._local_timezone>0: + header_content = header_content.replace(\ + '$TIMEZONE','UTC+'+str(config._local_timezone)) + elif config._local_timezone<0: + header_content = header_content.replace(\ + '$TIMEZONE','UTC'+str(config._local_timezone)) + + header_content = header_content.replace(\ + '$PROTOCOL_NUMBER',str(self.protocol_number)) + header_content = header_content.replace(\ + '$MODEL_NUMBER', str(self.model_number)) + header_content = header_content.replace(\ + '$FEATURE_NUMBER', str(self.feature_number)) + header_content = header_content.replace(\ + '$SERIAL_NUMBER', str(self.serial_number)) + + header_content = header_content.replace(\ + '$IXREADOUT', remove_linebreaks(self.ix_readout)) + header_content = header_content.replace(\ + '$RXREADOUT', remove_linebreaks(self.rx_readout)) + header_content = header_content.replace(\ + '$CXREADOUT', remove_linebreaks(self.cx_readout)) + + return(header_content) + + def format_content(self,timeutc_mean,timelocal_mean,temp_sensor,\ + freq_sensor,ticks_uC,sky_brightness): + # Format a string with data + date_time_utc_str = str(\ + timeutc_mean.strftime("%Y-%m-%dT%H:%M:%S"))+'.000' + date_time_local_str = str(\ + timelocal_mean.strftime("%Y-%m-%dT%H:%M:%S"))+'.000' + temp_sensor_str = str('%.2f' %temp_sensor) + ticks_uC_str = str('%.3f' %ticks_uC) + freq_sensor_str = str('%.3f' %freq_sensor) + sky_brightness_str = str('%.3f' %sky_brightness) + + formatted_data = \ + date_time_utc_str+";"+date_time_local_str+";"+temp_sensor_str+";"+\ + ticks_uC_str+";"+freq_sensor_str+";"+sky_brightness_str+"\n" + + return(formatted_data) + + def define_filenames(self): + # Filenames should follow a standard based on observatory name and date. + date_time_file = self.local_datetime(\ + self.read_datetime())-datetime.timedelta(hours=12) + date_file = date_time_file.date() + yearmonth = str(date_file)[0:7] + yearmonthday = str(date_file)[0:10] + + self.monthly_datafile = \ + config.monthly_data_directory+"/"+config._device_shorttype+\ + "_"+config._observatory_name+"_"+yearmonth+".dat" + #self.daily_datafile = \ + # config.daily_data_directory+"/"+config._device_shorttype+\ + # "_"+config._observatory_name+"_"+yearmonthday+".dat" + self.daily_datafile = \ + config.daily_data_directory+"/"+\ + yearmonthday.replace('-','')+'_120000_'+\ + config._device_shorttype+'-'+config._observatory_name+'.dat' + self.current_datafile = \ + config.current_data_directory+"/"+config._device_shorttype+\ + "_"+config._observatory_name+".dat" + + def save_data(self,formatted_data): + ''' + Save data to file and duplicate to current + data file (the one that will be ploted) + ''' + for each_file in [self.monthly_datafile,self.daily_datafile]: + if not os.path.exists(each_file): + datafile = open(each_file,'w') + datafile.write(self.standard_file_header()) + datafile.close() + + datafile = open(each_file,'a+') + datafile.write(formatted_data) + datafile.close() + + self.copy_file(self.daily_datafile,self.current_datafile) + + + def save_data_datacenter(self,formatted_data): + ''' + This function sends the data from this pysqm client to the central + node @ UCM. It saves the data there (only the SQM data file contents) + ''' + + # Connection details (hardcoded to avoid user changes) + DC_HOST = "muon.gae.ucm.es" + DC_PORT = 8739 + DEV_ID = str(config._device_id)+"_"+str(self.serial_number) + + def send_data(data): + try: + client = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + client.connect((DC_HOST, DC_PORT)) + client.sendall(data) + client.shutdown(socket.SHUT_RDWR) + client.close() + except: + return(0) + else: + return(1) + + def write_buffer(): + for data_line in self.DataBuffer[:]: + success = send_data(DEV_ID+";;D;;"+data_line) + if (success==1): self.DataBuffer.remove(data_line) + + return(success) + + ''' + Send the new file initialization to the datacenter + Appends the header to the buffer (it will be sent later) + ''' + + if (formatted_data=="NEWFILE"): + self.DataBuffer=[\ + hl+"\n" for hl in self.standard_file_header().split("\n")[:-1]] + + # Try to connect with the datacenter and send the header + success = send_data(DEV_ID+";;C;;") + success = write_buffer() + return(success) + else: + + ''' + Send the data to the datacenter + ''' + + # If the buffer is full, dont append more data. + if (len(self.DataBuffer)<10000): + self.DataBuffer.append(formatted_data) + + # Try to connect with the datacenter and send the data + success = write_buffer() + return(success) + + def save_data_mysql(self,formatted_data): + ''' + Use the Python MySQL API to save the + data to a database + ''' + mydb = None + values = formatted_data.split(';') + try: + ''' Start database connection ''' + mydb = _mysql.connect(\ + host = config._mysql_host, + user = config._mysql_user, + passwd = config._mysql_pass, + db = config._mysql_database, + port = config._mysql_port) + + ''' Insert the data ''' + mydb.query(\ + "INSERT INTO "+str(config._mysql_dbtable)+" VALUES (NULL,'"+\ + values[0]+"','"+values[1]+"',"+\ + values[2]+","+values[3]+","+\ + values[4]+","+values[5]+")") + except Exception as ex: + print((str(inspect.stack()[0][2:4][::-1])+\ + ' DB Error. Exception: %s' % str(ex))) + + if mydb != None: + mydb.close() + + def data_cache(self,formatted_data,number_measures=1,niter=0): + ''' + Append data to DataCache str. + If len(data)>number_measures, write to file + and flush the cache + ''' + try: + self.DataCache + except: + self.DataCache="" + + self.DataCache = self.DataCache+formatted_data + + if len(self.DataCache.split("\n"))>=number_measures+1: + self.save_data(self.DataCache) + self.DataCache = "" + print((str(niter)+'\t'+formatted_data[:-1])) + + def flush_cache(self): + ''' Flush the data cache ''' + self.save_data(self.DataCache) + self.DataCache = "" + + def copy_file(self,source,destination): + # Copy file content from source to dest. + fichero_source = open(source,'r') + contenido_source = fichero_source.read() + fichero_source.close() + # Create file and truncate it + fichero_destination = open(destination,'w') + fichero_destination.close() + # Write content + fichero_destination = open(destination,'r+') + fichero_destination.write(contenido_source) + fichero_destination.close() + + def remove_currentfile(self): + # Remove a file from the host + if os.path.exists(self.current_datafile): + os.remove(self.current_datafile) + + +class SQM(device): + def read_photometer(self,Nmeasures=1,PauseMeasures=2): + # Initialize values + temp_sensor = [] + flux_sensor = [] + freq_sensor = [] + ticks_uC = [] + Nremaining = Nmeasures + + # Promediate N measures to remove jitter + timeutc_initial = self.read_datetime() + while(Nremaining>0): + InitialDateTime = datetime.datetime.now() + + # Get the raw data from the photometer and process it. + raw_data = self.read_data(tries=10) + temp_sensor_i,freq_sensor_i,ticks_uC_i,sky_brightness_i = \ + self.data_process(raw_data) + + temp_sensor += [temp_sensor_i] + freq_sensor += [freq_sensor_i] + ticks_uC += [ticks_uC_i] + flux_sensor += [10**(-0.4*sky_brightness_i)] + Nremaining -= 1 + DeltaSeconds = (datetime.datetime.now()-InitialDateTime).total_seconds() + + # Just to show on screen that the program is alive and running + sys.stdout.write('.') + sys.stdout.flush() + + if (Nremaining>0): time.sleep(max(1,PauseMeasures-DeltaSeconds)) + + timeutc_final = self.read_datetime() + timeutc_delta = timeutc_final - timeutc_initial + + timeutc_mean = timeutc_initial+\ + datetime.timedelta(seconds=int(timeutc_delta.seconds/2.+0.5)) + timelocal_mean = self.local_datetime(timeutc_mean) + + # Calculate the mean of the data. + temp_sensor = filtered_mean(temp_sensor) + freq_sensor = filtered_mean(freq_sensor) + flux_sensor = filtered_mean(flux_sensor) + ticks_uC = filtered_mean(ticks_uC) + sky_brightness = -2.5*np.log10(flux_sensor) + + # Correct from offset (if cover is installed on the photometer) + #sky_brightness = sky_brightness+config._offset_calibration + + return(\ + timeutc_mean,timelocal_mean,\ + temp_sensor,freq_sensor,\ + ticks_uC,sky_brightness) + + def metadata_process(self,msg,sep=','): + # Separate the output array in items + msg = format_value(msg) + msg_array = msg.split(sep) + + # Get Photometer identification codes + self.protocol_number = int(format_value(msg_array[1])) + self.model_number = int(format_value(msg_array[2])) + self.feature_number = int(format_value(msg_array[3])) + self.serial_number = int(format_value(msg_array[4])) + + def data_process(self,msg,sep=','): + # Separate the output array in items + msg = format_value(msg) + msg_array = msg.split(sep) + + # Output definition characters + mag_char = 'm' + freq_char = 'Hz' + perc_char = 'c' + pers_char = 's' + temp_char = 'C' + + # Get the measures + sky_brightness = float(format_value(msg_array[1],mag_char)) + freq_sensor = float(format_value(msg_array[2],freq_char)) + ticks_uC = float(format_value(msg_array[3],perc_char)) + period_sensor = float(format_value(msg_array[4],pers_char)) + temp_sensor = float(format_value(msg_array[5],temp_char)) + + # For low frequencies, use the period instead + if freq_sensor<30 and period_sensor>0: + freq_sensor = 1./period_sensor + + return(temp_sensor,freq_sensor,ticks_uC,sky_brightness) + + def start_connection(self): + ''' Start photometer connection ''' + pass + + def close_connection(self): + ''' End photometer connection ''' + pass + + def reset_device(self): + ''' Restart connection''' + self.close_connection() + time.sleep(0.1) + #self.__init__() + self.start_connection() + + +class SQMLE(SQM): + def __init__(self): + ''' + Search the photometer in the network and + read its metadata + ''' + + try: + print(('Trying fixed device address %s ... ' %str(config._device_addr))) + self.addr = config._device_addr + self.port = 10001 + self.start_connection() + except: + print('Trying auto device address ...') + self.addr = self.search() + print(('Found address %s ... ' %str(self.addr))) + self.port = 10001 + self.start_connection() + + # Clearing buffer + print(('Clearing buffer ... |'), end=' ') + buffer_data = self.read_buffer() + print((buffer_data), end=' ') + print('| ... DONE') + print('Reading test data (ix,cx,rx)...') + time.sleep(1) + self.ix_readout = self.read_metadata(tries=10) + time.sleep(1) + self.cx_readout = self.read_calibration(tries=10) + time.sleep(1) + self.rx_readout = self.read_data(tries=10) + + def search(self): + ''' Search SQM LE in the LAN. Return its adress ''' + self.s = socket.socket(socket.AF_INET, socket.SOCK_DGRAM) + self.s.setblocking(False) + + if hasattr(socket,'SO_BROADCAST'): + self.s.setsockopt(socket.SOL_SOCKET, socket.SO_BROADCAST, 1) + + self.s.sendto("000000f6".decode("hex"), ("255.255.255.255", 30718)) + buf='' + starttime = time.time() + + print("Looking for replies; press Ctrl-C to stop.") + addr=[None,None] + while True: + try: + (buf, addr) = self.s.recvfrom(30) + if buf[3].encode("hex")=="f7": + print("Received from %s: MAC: %s" % \ + (addr, buf[24:30].encode("hex"))) + except: + #Timeout in seconds. Allow all devices time to respond + if time.time()-starttime > 3: + break + pass + + try: + assert(addr[0]!=None) + except: + print('ERR. Device not found!') + raise + else: + return(addr[0]) + + def start_connection(self): + ''' Start photometer connection ''' + self.s = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + self.s.settimeout(20) + self.s.connect((self.addr,int(self.port))) + #self.s.settimeout(1) + + def close_connection(self): + ''' End photometer connection ''' + self.s.setsockopt(\ + socket.SOL_SOCKET,\ + socket.SO_LINGER,\ + struct.pack('ii', 1, 0)) + + # Check until there is no answer from device + request = "" + r = True + while r: + r = self.read_buffer() + request += str(r) + self.s.close() + + def read_buffer(self): + ''' Read the data ''' + msg = None + try: msg = self.s.recv(256) + except: pass + return(msg) + + def reset_device(self): + ''' Connection reset ''' + #print('Trying to reset connection') + self.close_connection() + self.start_connection() + + def read_metadata(self,tries=1): + ''' Read the serial number, firmware version ''' + self.s.send('ix'.encode()) + time.sleep(1) + + read_err = False + msg = self.read_buffer().decode() + + # Check metadata + try: + # Sanity check + assert(len(msg)==_meta_len_ or _meta_len_==None) + assert("i," in msg) + self.metadata_process(msg) + except: + tries-=1 + read_err=True + + if (read_err==True and tries>0): + time.sleep(1) + self.reset_device() + time.sleep(1) + msg = self.read_metadata(tries) + if (msg!=-1): read_err=False + + # Check that msg contains data + if read_err==True: + print(('ERR. Reading the photometer!: %s' %str(msg))) + if (DEBUG): raise + return(-1) + else: + print(('Sensor info: '+str(msg)), end=' ') + return(msg) + + def read_calibration(self,tries=1): + ''' Read the calibration parameters ''' + self.s.send('cx'.encode()) + time.sleep(1) + + read_err = False + msg = self.read_buffer().decode() + + # Check caldata + try: + # Sanity check + assert(len(msg)==_cal_len_ or _cal_len_==None) + assert("c," in msg) + except: + tries-=1 + read_err=True + + if (read_err==True and tries>0): + time.sleep(1) + self.reset_device() + time.sleep(1) + msg = self.read_calibration(tries) + if (msg!=-1): read_err=False + + # Check that msg contains data + if read_err==True: + print(('ERR. Reading the photometer!: %s' %str(msg))) + if (DEBUG): raise + return(-1) + else: + print(('Calibration info: '+str(msg)), end=' ') + return(msg) + + def read_data(self,tries=1): + ''' Read the SQM and format the Temperature, Frequency and NSB measures ''' + self.s.send('rx'.encode()) + time.sleep(1) + + read_err = False + msg = self.read_buffer().decode() + + # Check data + try: + # Sanity check + assert(len(msg)==_data_len_ or _data_len_==None) + assert("r," in msg) + self.data_process(msg) + except: + tries-=1 + read_err=True + + if (read_err==True and tries>0): + time.sleep(1) + self.reset_device() + time.sleep(1) + msg = self.read_data(tries) + if (msg!=-1): read_err=False + + # Check that msg contains data + if read_err==True: + print(('ERR. Reading the photometer!: %s' %str(msg))) + if (DEBUG): raise + return(-1) + else: + if (DEBUG): print(('Data msg: '+str(msg))) + return(msg) + + +class SQMLU(SQM): + def __init__(self): + ''' + Search the photometer and + read its metadata + ''' + + try: + print(('Trying fixed device address %s ... ' %str(config._device_addr))) + self.addr = config._device_addr + self.bauds = 115200 + self.start_connection() + except: + print('Trying auto device address ...') + self.addr = self.search() + print(('Found address %s ... ' %str(self.addr))) + self.bauds = 115200 + self.start_connection() + + # Clearing buffer + print(('Clearing buffer ... |'), end=' ') + buffer_data = self.read_buffer() + print((buffer_data), end=' ') + print('| ... DONE') + print('Reading test data (ix,cx,rx)...') + time.sleep(1) + self.ix_readout = self.read_metadata(tries=10) + time.sleep(1) + self.cx_readout = self.read_calibration(tries=10) + time.sleep(1) + self.rx_readout = self.read_data(tries=10) + + + def search(self): + ''' + Photometer search. + Name of the port depends on the platform. + ''' + ports_unix = ['/dev/ttyUSB'+str(num) for num in range(100)] + ports_win = ['COM'+str(num) for num in range(100)] + + os_in_use = sys.platform + + if os_in_use == 'linux2': + print('Detected Linux platform') + ports = ports_unix + elif os_in_use == 'win32': + print('Detected Windows platform') + ports = ports_win + + used_port = None + for port in ports: + conn_test = serial.Serial(port, 115200, timeout=1) + conn_test.write('ix'.encode) + if conn_test.readline()[0] == 'i': + used_port = port + break + + try: + assert(used_port!=None) + except: + print('ERR. Device not found!') + raise + else: + return(used_port) + + def start_connection(self): + '''Start photometer connection ''' + + self.s = serial.Serial(self.addr, 115200, timeout=2) + + def close_connection(self): + ''' End photometer connection ''' + # Check until there is no answer from device + request = "" + r = True + while r: + r = self.read_buffer() + request += str(r) + + self.s.close() + + def reset_device(self): + ''' Connection reset ''' + #print('Trying to reset connection') + self.close_connection() + self.start_connection() + + def read_buffer(self): + ''' Read the data ''' + msg = None + try: msg = self.s.readline() + except: pass + return(msg) + + def read_metadata(self,tries=1): + ''' Read the serial number, firmware version ''' + self.s.write("ix".encode()) + time.sleep(1) + + read_err = False + msg = self.read_buffer().decode() + + # Check metadata + try: + # Sanity check + assert(len(msg)==_meta_len_ or _meta_len_==None) + assert("i," in msg) + self.metadata_process(msg) + except: + tries-=1 + read_err=True + + if (read_err==True and tries>0): + time.sleep(1) + self.reset_device() + time.sleep(1) + msg = self.read_metadata(tries) + if (msg!=-1): read_err=False + + # Check that msg contains data + if read_err==True: + print(('ERR. Reading the photometer!: %s' %str(msg))) + if (DEBUG): raise + return(-1) + else: + print(('Sensor info: '+str(msg)), end=' ') + return(msg) + + def read_calibration(self,tries=1): + ''' Read the calibration data ''' + self.s.write('cx'.encode()) + time.sleep(1) + + read_err = False + msg = self.read_buffer().decode() + + # Check caldata + try: + # Sanity check + assert(len(msg)==_cal_len_ or _cal_len_==None) + assert("c," in msg) + except: + tries-=1 + read_err=True + + if (read_err==True and tries>0): + time.sleep(1) + self.reset_device() + time.sleep(1) + msg = self.read_calibration(tries) + if (msg!=-1): read_err=False + + # Check that msg contains data + if read_err==True: + print(('ERR. Reading the photometer!: %s' %str(msg))) + if (DEBUG): raise + return(-1) + else: + print(('Calibration info: '+str(msg)), end=' ') + return(msg) + + def read_data(self,tries=1): + ''' Read the SQM and format the Temperature, Frequency and NSB measures ''' + self.s.write('rx'.encode()) + time.sleep(1) + + read_err = False + msg = self.read_buffer().decode() + + # Check data + try: + # Sanity check + assert(len(msg)==_data_len_ or _data_len_==None) + assert("r," in msg) + self.data_process(msg) + except: + tries-=1 + read_err=True + + if (read_err==True and tries>0): + time.sleep(1) + self.reset_device() + time.sleep(1) + msg = self.read_data(tries) + if (msg!=-1): read_err=False + + # Check that msg contains data + if read_err==True: + print(('ERR. Reading the photometer!: %s' %str(msg))) + if (DEBUG): raise + return(-1) + else: + if (DEBUG): print(('Data msg: '+str(msg))) + return(msg) + diff --git a/pysqm/settings.py b/pysqm/settings.py new file mode 100644 index 0000000000000000000000000000000000000000..00e995892a19230aa0357f16b5375cc9483c6a41 --- /dev/null +++ b/pysqm/settings.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python + +''' +PySQM plotting program +____________________________ + +Copyright (c) Mireia Nievas <mnievas[at]ucm[dot]es> + +This file is part of PySQM. + +PySQM is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +PySQM is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with PySQM. If not, see <http://www.gnu.org/licenses/>. +____________________________ +''' + +import os,sys + + +class ArgParser: + def __init__(self,inputfile=False): + self.parse_arguments(inputfile) + + def parse_arguments(self,inputfile): + import argparse + # Return config filename + self.parser = argparse.ArgumentParser() + self.parser.add_argument('-c', '--config', default="config.py") + if (inputfile): + self.parser.add_argument('-i', '--input', default=None) + args = self.parser.parse_args() + vars(self).update(args.__dict__) + + def print_help(self): + self.parser.print_help() + + +class ConfigFile: + def __init__(self, path="config.py"): + # Guess the selected dir and config filename + # Should accept: + # - absolute path (inc. filename) + # - relative path (inc. filename) + # - absolute path (exc. filename) + # - relative path (exc. filename) + # - shortcouts like ~ . etc + self.path = path + self.config = None + + def read_config_file(self,path): + # Get the absolute path + abspath = os.path.abspath(path) + # Is a dir? Then add config.py (default filename) + if os.path.isdir(abspath): + abspath += "/config.py" + # split directory and filename + directory = os.path.dirname(abspath) + filename = os.path.basename(abspath) + + old_syspath = sys.path + sys.path.append(directory) + import config + self.config = config + +# Create an object (by default empty) accessible from everywhere +# After read_config_file is called, GlobalConfig.config will be accessible +GlobalConfig = ConfigFile() diff --git a/setup.py b/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..75a6ec9907e283ef1c9bcacc136964601b90c6fe --- /dev/null +++ b/setup.py @@ -0,0 +1,27 @@ +#!/usr/bin/env python + +from distutils.core import setup + +setup(name='pysqm', + version='3.1', + maintainer='Thorsten Alteholz', + maintainer_email='python@alteholz.de', + url='https://github.com/alteholz/PySQM', + license='GPLv3', + description='SQM reading and plotting software', + packages=['pysqm'], + install_requires=['pyephem','numpy','matplotlib'], + classifiers=[ + "Programming Language :: C", + "Programming Language :: Cython", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: Implementation :: CPython", + 'Development Status :: 3 - Alpha', + "Environment :: Other Environment", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: GNU General Public License (GPL)", + "Operating System :: OS Independent", + "Topic :: Scientific/Engineering :: Astronomy", + ], + long_description=open('README.txt').read() + )