source.py

from datetime import datetime
import csv
from collections import OrderedDict
import math
from statistics import median
from statistics import mean
from matplotlib import pyplot as plt
import numpy as np

# DATE formats
geiger_format = '%Y-%m-%d %H:%M:%S'
alt_format    = '%Y-%m-%dT%H:%M:%S.%fZ'

# Read the data files

geiger_data = OrderedDict()
with open('geiger_log.csv') as csv_file:
    csv_reader = csv.reader(csv_file,delimiter=',')
    linecount = 0
    for row in csv_reader:
        if linecount > 0:
            #print(f'{row[0]}   {row[6]}')
            date = datetime.strptime(row[0],geiger_format).timestamp()
            geiger_data[date] = row[6]
        linecount += 1

alt_data = OrderedDict()
with open('gps_log.csv') as csv_file:
    csv_reader = csv.reader(csv_file,delimiter=',')
    linecount = 0
    for row in csv_reader:
        if linecount>0:
            date = datetime.strptime(row[0],alt_format).timestamp()
            alt_data[date] = row[3]
        linecount += 1


# Create joined data
joined_data_geiger = OrderedDict()
joined_data_alt = OrderedDict()
start_time = datetime.strptime('2019-09-14 12:20:00',geiger_format).timestamp()
end_time = datetime.strptime('2019-09-14 17:20:00',geiger_format).timestamp()
for row in geiger_data:
    if row > start_time and row < end_time:
        # Find the closest smaller and higher altitude times, average them
        before_alt_time = 0
        after_alt_time = 0
        keyList=sorted(alt_data.keys())
        for i,v in enumerate(keyList):
            if v>=row:
                scale = (row-keyList[i-1])/(keyList[i]-keyList[i-1])
                value1 = float(alt_data[keyList[i-1]])
                value2 = float(alt_data[keyList[i]])
                diff = value2-value1
                diff *= scale
                interp_value = value1+scale
                time = datetime.fromtimestamp(row-60).strftime("%H:%M")
                joined_data_geiger[time]=float(geiger_data[row])
                joined_data_alt[time]=interp_value
                break

# Subroutine for calculating the average line in the graph
def averages(x_array,y_array,min,max,binsize):
    """Calculate the average of the binned data. The bins are along the x axis"""
    bincount = 1+math.floor((max-min)/binsize)
    total = np.zeros(bincount)
    count = np.zeros(bincount)
    for index in range(len(x_array)):
        y = y_array[index]
        x = x_array[index]
        bin = math.floor((x-min)/binsize)
        total[bin] = total[bin]+y
        count[bin] = count[bin]+1

    for i in range(bincount):
        if count[i]>0:
            total[i]=total[i]/count[i]

    return total

# Plot the time log graph
from matplotlib import pyplot as plt
%matplotlib inline

font = {'family' : 'normal',
        'weight' : 'normal',
        'size'   : 18}

plt.rc('font', **font)
fig=plt.figure()

color = 'tab:blue'
ax=fig.add_axes([0,0,2,2])
ax.set_ylabel('altitude, m',color=color)
ax.plot(list(joined_data_alt.keys()),list(joined_data_alt.values()),'-',color=color)
ax.tick_params(axis='y', labelcolor=color)

ax2 = ax.twinx()  # instantiate a second axes that shares the same x-axis

color = 'tab:red'
ax2.set_ylabel('dose rate, μSv/h', color=color)  # we already handled the x-label with ax1
ax2.plot(list(joined_data_geiger.keys()),list(joined_data_geiger.values()),color=color)
ax2.tick_params(axis='y', labelcolor=color)


ax.set_title('Altitude and radiation dose rate during two fligths from FLR via FRA to GOT ')
ax.set_xlabel('time (UTC)')
every_nth = 10
for n, label in enumerate(ax.xaxis.get_ticklabels()):
    if (n) % every_nth != 0:
        label.set_visible(False)
plt.show()

# Plot the figure showing radiation vs. altitude
plt.rc('font', **font)
fig=plt.figure()

x = list(joined_data_alt.values())
y = list(joined_data_geiger.values())
binsize =  1000

color = 'tab:red'
ax=fig.add_axes([0,0,1,1])
ax.set_ylabel('μSv/h',color=color)
ax.plot(x,y,'.',color=color)
ax.tick_params(axis='y', labelcolor=color)
xmin = math.floor(min(x))
xmax = math.floor(max(x))
mean = averages(x,y,xmin,xmax,binsize)
ax.plot( list(range(int(xmin+binsize/2),int(xmax+1+binsize/2),binsize)),mean,'-',color='black')


ax.set_title('Radiation dose rate versus altitude')
ax.set_xlabel('altitude, m')
plt.show()