#---------------------------ooo0ooo---------------------------
#       EFM monitor v5.1
#       updated 17/04/2021
#		Ian Robinson
#       www.starfishprime.co.uk
#---------------------------Notes---------------------------
#
#   each data point     = 37 bytes
#   at 1 min intervals  = 20Mb p.a
#   dataFile format (txt) = datetime .. B Field nT .. temperature celcius
#
#   older data is periodically stripped out and archived
#       to save memory and speed processing
#
#   ActiveData.dat[0][i] = date
#   ActiveData.dat[1][i] = field
#   ActiveData.dat[2][i] = temperature
#
#---------------------------ooo0ooo---------------------------

import serial
import time

import string
import matplotlib
matplotlib.use('Agg') #prevent use of Xwindows
from matplotlib import pyplot as plt
from threading import Thread
import statistics   #used by median filter
import os
import gc

import numpy as np
import time
import array
from array import array
from scipy import signal
from obspy import UTCDateTime, read, Trace, Stream
from shutil import copyfile



#---------------------------ooo0ooo---------------------------
def createMSeed(DataArray, StartDateTime, EndDateTime, nDataPoints, stationid, stationchannel, location):


    trueSampleFreq = float(nDataPoints) / (EndDateTime - StartDateTime)

    # Fill header attributes
    stats = {'network': 'IR', 'station': stationid, 'location': location,
         'channel': stationchannel, 'npts': nDataPoints, 'sampling_rate': trueSampleFreq,
         'mseed': {'dataquality': 'D'}}
    # set current time
    stats['starttime'] = StartDateTime
    st = Stream([Trace(data=DataArray[0:nDataPoints], header=stats)])
    return st
#---------------------------ooo0ooo---------------------------
def SaveDataMSeed(st, StartDateTime, stationid, stationchannel, location):


  Year = str(StartDateTime.year)
  Month = str(StartDateTime.month)
  Day = str(StartDateTime.day)
  Hour = str(StartDateTime.hour)
  Minute = str(StartDateTime.minute)

  yearDir =  'Data' + '/' + Year

  FileName = str(Year)+'-'+str(Month)+'-'+str(Day)+'__'+ Hour+':'+Minute +'__'+str(stationid)+'-'+str(stationchannel)+'-'+str(location)+'.mseed'
  
  here = os.path.dirname(os.path.realpath(__file__))

  try:
    os.makedirs(yearDir)
  except OSError:
    if not os.path.isdir(yearDir):
        raise

  FilePath = os.path.join(here, yearDir, FileName)
  
  dataFile = open(FilePath, 'wb')
  st.write(dataFile,format='MSEED',  encoding=4, reclen=256)
  dataFile.close()
  print ('Data write of active data completed')

#---------------------------ooo0ooo---------------------------
def SaveDailyHSDMSeed(st, StartDateTime, stationid, stationchannel, location):


  Year = str(StartDateTime.year)
  Month = str(StartDateTime.month)
  Day = str(StartDateTime.day)
  Hour = str(StartDateTime.hour)
  Minute = str(StartDateTime.minute)

  yearDir =  'Data' + '/' + Year

  FileName = str(Year)+'-'+str(Month)+'-'+str(Day)+'__'+ Hour+':'+Minute +'__'+str(stationid)+'-'+str(stationchannel)+'-'+str(location)+'_daily_HSD.mseed'
  
  here = os.path.dirname(os.path.realpath(__file__))

  try:
    os.makedirs(yearDir)
  except OSError:
    if not os.path.isdir(yearDir):
        raise

  FilePath = os.path.join(here, yearDir, FileName)
  
  dataFile = open(FilePath, 'wb')
  st.write(dataFile,format='MSEED',  encoding=4, reclen=256)
  dataFile.close()
  print ('Data write of active data completed')


#---------------------------ooo0ooo--------------------------  
	
def SaveDataArray(dataArray, f_name):
	#simple routine to save an array as a text file - used for testing purposes
	
	file = open(f_name, "w+")
	for i in dataArray:
		file.write(str(i) +'\n')

	file.close()
	return
#---------------------------ooo0ooo---------------------------
def plotTodaysBField(DataArray, StartDateTime, EndDateTime, nDataPoints):

  print('plotting DayPlot',nDataPoints)
  
  if (nDataPoints > 20):
    ydata=DataArray[3:nDataPoints]
    

    data=removeGlitches(ydata)
    data=medianFilter(ydata, len(ydata))

    zeroPoint = np.mean(data)
    data=data-zeroPoint #(de-mean the data)
    
    
    yearNow = StartDateTime.year
    monthNow = StartDateTime.month
    dayNow = StartDateTime.day
    
    dateString = str(yearNow) + '-' + str(monthNow) + '-' + str(dayNow)
    filename1 = ('Plots/Today.svg')
    filename2 = 'Plots/' + dateString + '.svg'
    


    YAxisRange=np.amax(data)*1.1
    
    startHour=StartDateTime.hour
    startMinute=StartDateTime.minute
    timeSampled = EndDateTime-StartDateTime  #length of data in seconds
    graphStart = float(startHour) + float(startMinute/60)
    graphEnd = graphStart + (timeSampled/3600.0)
    xValues= np.linspace(graphStart, graphEnd, len(ydata))
    upDated = str(UTCDateTime())
        
    fig = plt.figure(figsize=(12,4))
    xAxisText="Time  updated - "+ upDated +" UTC"
    plt.title('Magnetic Field E-W, Guisborough - '+ dateString)
    plt.xlabel(xAxisText)
    zeroLabel=str('%0.5g' % (zeroPoint))
    plt.ylabel('$\Delta$Flux Density - nT     0.0='+zeroLabel+'nT')
    plt.plot(xValues, data,  marker='None',    color = 'darkolivegreen')  
    plt.xlim(0, 24.1)
    plt.xticks(np.arange(0, 24.1, 3.0))
    plt.grid(True)
    plt.savefig(filename1)
    
    copyfile('Plots/Today.svg', filename2)   
    plt.close('all')     
    

  return
  
#---------------------------ooo0ooo---------------------------
def plotTodaysHSD(dailyStDevs,nDailyHSD,StartDateTime,EndDateTime,prevHrCollected):
	##on the first day of running the hsd graph is 1 hour behind the normal plot - this needs sorting
	print('plotting hsdPlot',nDailyHSD)

	try:
		if (prevHrCollected == True):
			yearNow = StartDateTime.year
			monthNow = StartDateTime.month
			dayNow = StartDateTime.day
			dateString = str(yearNow) + '-' + str(monthNow) + '-' + str(dayNow)
			filename1 = ('Plots/TodayHSD.svg')
			filename2 = ('Plots/' + dateString + '_HSD.svg')
			YAxisRange=(np.amax(dailyStDevs))*1.5

			startHour=StartDateTime.hour
			startMinute=StartDateTime.minute
			timeSampled = EndDateTime-StartDateTime  #length of data in seconds
			graphStart = float(startHour) + float(startMinute/60)
			graphEnd = graphStart + (timeSampled/3600.0)

			upDated = str(UTCDateTime())

			yValues=dailyStDevs[0:(nDailyHSD-1)]		
			xValues= np.linspace(graphStart, graphEnd, len(yValues))
			
			
			fig = plt.figure(figsize=(12,4))
			xAxisText=("Time  updated - "+ upDated +" UTC")
			plt.title('Hourly Standard Deviations - '+ dateString)
			plt.xlabel(xAxisText)
		
			plt.ylabel('$\sigma$ - nT')

			plt.plot(xValues, yValues,  marker='None',    color = 'darkolivegreen')  
			plt.xlim(0, 24.1)
			plt.ylim(0, YAxisRange)
			plt.xticks(np.arange(0, 24.1, 3.0))
			plt.grid(True)
			
			plt.savefig(filename1)
			copyfile(filename1, filename2) 
			plt.close('all') 

	except ValueError:
		print ('failed to plot weeklyhsd graph')	
    
	
	
	return
#---------------------------ooo0ooo---------------------------
def plotWeeksHSD(weeklyStDevs,nWeeklyHSD,weekStartDateTime,EndDateTime,prevHrCollected):

	try:
		if (prevHrCollected == True):
			yearNow = weekStartDateTime.year
			monthNow = weekStartDateTime.month
			dayNow = weekStartDateTime.day
			dateString = str(yearNow) + '-' + str(monthNow) + '-' + str(dayNow)
			filename1 = ('Plots/ThisWeekHSD.svg')
			filename2 = ('Plots/' + dateString + '_weekly_HSD.svg')
			YAxisRange=(np.amax(weeklyStDevs))*1.5

			startDay=weekStartDateTime.isoweekday()
			startHour=weekStartDateTime.hour
			startMinute=weekStartDateTime.minute
			timeSampled = EndDateTime-weekStartDateTime  #length of data in seconds
			graphStart = float(startDay-1) + (float(startHour)/24) + float(startMinute/1440)
			graphEnd = graphStart + (timeSampled/86400.0)


			upDated = str(UTCDateTime())
			dateString = str(yearNow) + '-' + str(monthNow) + '-' + str(dayNow)
			yValues=weeklyStDevs[0:(nWeeklyHSD-1)]		
			xValues= np.linspace(graphStart, graphEnd, len(yValues))
			
			
			fig = plt.figure(figsize=(12,4))
			xAxisText=("Time  updated - "+ upDated +" UTC")
			plt.title('Earth Mag Field Weekly Standard Deviations - wc- '+ dateString)
			plt.xlabel(xAxisText)
		
			plt.ylabel('$\sigma$ - nT')

			plt.plot(xValues, yValues,  marker='None',    color = 'darkolivegreen')  
			plt.xlim(0, 7.01)
			plt.ylim(0, YAxisRange)
			#plt.xticks(np.arange(0, 7.01, 1.0))
			#ticks = list(range(0, 7, 1))
			ticks=np.arange(0, 7.01, 1.0)
			labels = "Mon Tues Weds Thurs Fri Sat Sun".split()
			plt.xticks(ticks, labels)
			plt.grid(True)
			
			plt.savefig(filename1)
			copyfile(filename1, filename2)  
			plt.close('all') 

	except ValueError:
		print ('failed to plot weeklyhsd graph')	
    
	
	
	return

#---------------------------ooo0ooo---------------------------
def medianFilter(ydata, z):

  dataNew=ydata

  for i in range (2,z-2):
    medData=ydata[i]
    dataChunk=[ydata[i-2],ydata[i-1],ydata[i+1],ydata[i+2]]
    medianDatum=statistics.median(dataChunk)
    dataNew[i]=medianDatum
  return dataNew

#---------------------------ooo0ooo---------------------------
def removeGlitches(ydata):
  #to maintain continuity of data-line miss-reads are assigned B=0.0
  z=len(ydata)
  dataNew=ydata

  for i in range (0,z-4):
    if (abs(ydata[i]) < 0.1):
        dataChunk=[ydata[i],ydata[i+1],ydata[i+2],ydata[i+3]]
        dataNew[i]=statistics.median(dataChunk)
  
  dataNew[(z-3):(z-1)]=ydata[(z-3):(z-1)]
  return dataNew
#---------------------------ooo0ooo---------------------------
def runningMeanFast(x, N):
  return np.convolve(x, np.ones((N,))/N)[(N-1):]
#---------------------------ooo0ooo---------------------------
def SaveAndPlot(DataArray, StartDateTime, EndDateTime, nDataPoints, stationid,stationchannel, location, dailyStDevs,nDailyHSD,prevHrCollected,weeklyStDevs,nWeeklyHSD,weekStartDateTime,):
	st = createMSeed(DataArray, StartDateTime, EndDateTime, nDataPoints, stationid, stationchannel, location)
	SaveDataMSeed(st, StartDateTime, stationid, stationchannel, location)
	st = createMSeed(dailyStDevs, StartDateTime, EndDateTime, nDailyHSD, stationid, stationchannel, location)
	SaveDailyHSDMSeed(st, StartDateTime, stationid, stationchannel, location)
	plotTodaysBField(DataArray, StartDateTime, EndDateTime, nDataPoints)
	plotTodaysHSD(dailyStDevs,nDailyHSD,StartDateTime,EndDateTime,prevHrCollected)
	plotWeeksHSD(weeklyStDevs,nWeeklyHSD,weekStartDateTime,EndDateTime,prevHrCollected)
	return
#---------------------------ooo0ooo---------------------------
#---------------------------ooo0ooo---------------------------

#-------------------------------------------------------------------
#---------------------------ooo0ooo---------------------------
#                                                            #
#                   main body                                #
#                                                            #
#---------------------------ooo0ooo--------------------------#

###*******FGM-3 Sensor Values**********
mSens = 0.118
cSens = 14.353
###************************************
###************************************

###*******TM36 temp probe conversion  V to degrees C*********
tempScaler=0.0100
tempOffset = 0.5000
###************************************
#-- station parameters
stationid = 'STARF'
# channel L=low Freq, F=magnetometer E = east/west orientation
stationchannel = 'LFE'  #see SEED format documentation
location = '01'  # 2 digit code to identify specific sensor rig


#-------------create folders if absent------------------------
os.chdir('/home/pi/EFM')
try:
    os.makedirs('Plots/')
except OSError:
    if not os.path.isdir('Plots/'):
        raise
###---------------------------ooo0ooo---------------------------




DataFileLengthSeconds = 86400

# time between samples in seconds      
SamplingPeriod = 2.00


# time between datasaves and plots in seconds, 
DataWriteInterval = 1000

# time (seconds) between determination of Hourly Standard Deviations 
hsdInterval=60

# -- create numpy array to store dailydata
TargetNoSamples = int((DataFileLengthSeconds*1.1)/SamplingPeriod)
DataArray=np.zeros([TargetNoSamples], dtype=np.float32)
nDataPoints = 0
FirstDataReading=True

# -- create numpy array for previous 1 hour data for determination
# -- of HSD (Hourly Standard Deviation
TargetHourlyHSDSamples = int((3600)/SamplingPeriod)
prevHrData=np.zeros([TargetHourlyHSDSamples], dtype=np.float32)
headPointer=0  #implementing prevHrData as a circular queue
prevHrCollected = False  # true when a full hour's data has been collected

# -- create Numpy array for day's hourly standard deviations  
TargetDailyHSDSamples = int((86400*1.05)/hsdInterval)
dailyStDevs=np.zeros([TargetDailyHSDSamples], dtype=np.float32)
nDailyHSD=0

# -- create Numpy array to store week's hourly standard deviations  
TargetWeeklyHSDSamples = int((608400*1.05)/hsdInterval)
weeklyStDevs=np.zeros([TargetWeeklyHSDSamples], dtype=np.float32)
nWeeklyHSD=0


time.sleep(30)  #give Arduino time to boot
ser = serial.Serial('/dev/ttyACM0', 57600)
ser.flushInput()
time.sleep(2)

weekStartDateTime=UTCDateTime()
StartDateTime=weekStartDateTime
lastSaveTime=StartDateTime
lastHSDcalcTime=StartDateTime



while 1:
	time.sleep(SamplingPeriod)
	
	try:
		freq = np.float32(ser.readline())   # converts the freq number string to a floating point number
		if (freq > 10000):
			Bfield = (((1000000.0/freq)-cSens)/mSens)*1000.0# convert freq to field in nTeslas
		else:
			Bfield=0.0

		if (FirstDataReading == False):  #reject first reading which is often spurious
			DataArray[nDataPoints]=Bfield
			nDataPoints = nDataPoints +1
			prevHrData[headPointer]=Bfield
			headPointer=headPointer+1
      
			if (headPointer >= (TargetHourlyHSDSamples-1)):
				prevHrCollected = True #can now start to calculate s.dev
				headPointer=0

		else:
			FirstDataReading=False # first reading rejected, subsequent will be accepted
			StartDateTime=UTCDateTime()
			lastSaveTime = UTCDateTime()

	except ValueError:
		print ('failed to read serial line correctly --trying again')



      
	if (StartDateTime.day != UTCDateTime().day): # entered a new day
		EndDateTime=UTCDateTime()
		threadSaveAndPlot = Thread(target=SaveAndPlot, args=(DataArray, StartDateTime, EndDateTime, nDataPoints, stationid,stationchannel, location, dailyStDevs,nDailyHSD,prevHrCollected,weeklyStDevs,nWeeklyHSD,weekStartDateTime,))
		threadSaveAndPlot.start()
		
		if weekStartDateTime.isoweekday()==1: #Monday-- new week starts
			#clear theweekly data array
			weeklyStDevs=np.zeros([TargetWeeklyHSDSamples], dtype=np.float32)
			nWeeklyHSD=0
			
		DataArray=np.zeros(TargetNoSamples, np.float32)  # zero data array for new day
		dailyStDevs=np.zeros([TargetDailyHSDSamples], dtype=np.float32)
		nDailyHSD=0
		StartDateTime=UTCDateTime()
		lastSaveTime=StartDateTime
		nDataPoints = 0
		gc.collect

	if ((UTCDateTime() - lastSaveTime) > DataWriteInterval):
		EndDateTime=UTCDateTime()
		threadSaveAndPlot = Thread(target=SaveAndPlot, args=(DataArray, StartDateTime, EndDateTime, nDataPoints, stationid,stationchannel, location, dailyStDevs,nDailyHSD,prevHrCollected,weeklyStDevs,nWeeklyHSD,weekStartDateTime,))
		threadSaveAndPlot.start()
		lastSaveTime = UTCDateTime()
		
	if((UTCDateTime() - lastHSDcalcTime) > hsdInterval):
		#perform numpy calc of standard deviation of previous hours data
		if (prevHrCollected == True):
			prevHourHSD=np.std(prevHrData[0:(TargetHourlyHSDSamples-1)])
			dailyStDevs[nDailyHSD]= prevHourHSD
			weeklyStDevs[nWeeklyHSD]= prevHourHSD
			nDailyHSD=nDailyHSD+1
			nWeeklyHSD=nWeeklyHSD+1
		lastHSDcalcTime= UTCDateTime()
		


      

