SHARE
TWEET

source.py

a guest Sep 18th, 2019 122 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. from datetime import datetime
  2. import csv
  3. from collections import OrderedDict
  4. import math
  5. from statistics import median
  6. from statistics import mean
  7. from matplotlib import pyplot as plt
  8. import numpy as np
  9.  
  10. # DATE formats
  11. geiger_format = '%Y-%m-%d %H:%M:%S'
  12. alt_format    = '%Y-%m-%dT%H:%M:%S.%fZ'
  13.  
  14. # Read the data files
  15.  
  16. geiger_data = OrderedDict()
  17. with open('geiger_log.csv') as csv_file:
  18.     csv_reader = csv.reader(csv_file,delimiter=',')
  19.     linecount = 0
  20.     for row in csv_reader:
  21.         if linecount > 0:
  22.             #print(f'{row[0]}   {row[6]}')
  23.             date = datetime.strptime(row[0],geiger_format).timestamp()
  24.             geiger_data[date] = row[6]
  25.         linecount += 1
  26.  
  27. alt_data = OrderedDict()
  28. with open('gps_log.csv') as csv_file:
  29.     csv_reader = csv.reader(csv_file,delimiter=',')
  30.     linecount = 0
  31.     for row in csv_reader:
  32.         if linecount>0:
  33.             date = datetime.strptime(row[0],alt_format).timestamp()
  34.             alt_data[date] = row[3]
  35.         linecount += 1
  36.    
  37.  
  38. # Create joined data
  39. joined_data_geiger = OrderedDict()
  40. joined_data_alt = OrderedDict()
  41. start_time = datetime.strptime('2019-09-14 12:20:00',geiger_format).timestamp()
  42. end_time = datetime.strptime('2019-09-14 17:20:00',geiger_format).timestamp()
  43. for row in geiger_data:
  44.     if row > start_time and row < end_time:
  45.         # Find the closest smaller and higher altitude times, average them
  46.         before_alt_time = 0
  47.         after_alt_time = 0
  48.         keyList=sorted(alt_data.keys())
  49.         for i,v in enumerate(keyList):
  50.             if v>=row:
  51.                 scale = (row-keyList[i-1])/(keyList[i]-keyList[i-1])
  52.                 value1 = float(alt_data[keyList[i-1]])
  53.                 value2 = float(alt_data[keyList[i]])
  54.                 diff = value2-value1
  55.                 diff *= scale
  56.                 interp_value = value1+scale
  57.                 time = datetime.fromtimestamp(row-60).strftime("%H:%M")
  58.                 joined_data_geiger[time]=float(geiger_data[row])
  59.                 joined_data_alt[time]=interp_value
  60.                 break
  61.  
  62. # Subroutine for calculating the average line in the graph
  63. def averages(x_array,y_array,min,max,binsize):
  64.     """Calculate the average of the binned data. The bins are along the x axis"""
  65.     bincount = 1+math.floor((max-min)/binsize)
  66.     total = np.zeros(bincount)
  67.     count = np.zeros(bincount)
  68.     for index in range(len(x_array)):
  69.         y = y_array[index]
  70.         x = x_array[index]
  71.         bin = math.floor((x-min)/binsize)
  72.         total[bin] = total[bin]+y
  73.         count[bin] = count[bin]+1
  74.    
  75.     for i in range(bincount):
  76.         if count[i]>0:
  77.             total[i]=total[i]/count[i]
  78.        
  79.     return total
  80.  
  81. # Plot the time log graph
  82. from matplotlib import pyplot as plt
  83. %matplotlib inline
  84.  
  85. font = {'family' : 'normal',
  86.         'weight' : 'normal',
  87.         'size'   : 18}
  88.  
  89. plt.rc('font', **font)
  90. fig=plt.figure()
  91.  
  92. color = 'tab:blue'
  93. ax=fig.add_axes([0,0,2,2])
  94. ax.set_ylabel('altitude, m',color=color)
  95. ax.plot(list(joined_data_alt.keys()),list(joined_data_alt.values()),'-',color=color)
  96. ax.tick_params(axis='y', labelcolor=color)
  97.  
  98. ax2 = ax.twinx()  # instantiate a second axes that shares the same x-axis
  99.  
  100. color = 'tab:red'
  101. ax2.set_ylabel('dose rate, μSv/h', color=color)  # we already handled the x-label with ax1
  102. ax2.plot(list(joined_data_geiger.keys()),list(joined_data_geiger.values()),color=color)
  103. ax2.tick_params(axis='y', labelcolor=color)
  104.  
  105.  
  106. ax.set_title('Altitude and radiation dose rate during two fligths from FLR via FRA to GOT ')
  107. ax.set_xlabel('time (UTC)')
  108. every_nth = 10
  109. for n, label in enumerate(ax.xaxis.get_ticklabels()):
  110.     if (n) % every_nth != 0:
  111.         label.set_visible(False)
  112. plt.show()
  113.  
  114. # Plot the figure showing radiation vs. altitude
  115. plt.rc('font', **font)
  116. fig=plt.figure()
  117.  
  118. x = list(joined_data_alt.values())
  119. y = list(joined_data_geiger.values())
  120. binsize =  1000
  121.  
  122. color = 'tab:red'
  123. ax=fig.add_axes([0,0,1,1])
  124. ax.set_ylabel('μSv/h',color=color)
  125. ax.plot(x,y,'.',color=color)
  126. ax.tick_params(axis='y', labelcolor=color)
  127. xmin = math.floor(min(x))
  128. xmax = math.floor(max(x))
  129. mean = averages(x,y,xmin,xmax,binsize)
  130. ax.plot( list(range(int(xmin+binsize/2),int(xmax+1+binsize/2),binsize)),mean,'-',color='black')
  131.  
  132.  
  133. ax.set_title('Radiation dose rate versus altitude')
  134. ax.set_xlabel('altitude, m')
  135. plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top