Radiocarbon emissions influence on background ratio 14C to 12C
- Shift in ration between radiocarbon and carbon because of nuclear power plants and fuel reprocessing stations
- Access and display radiocarbon emission year 2016 based on data from the supplemental material of the paper Zazzeri et al. 2018
- Influence of nuclear power plants and fuel reprocessing stations: Timeseries by station.
- Influence of nuclear power plants and fuel reprocessing stations: Maps and bar graph.
Background
There are several different isotopes of carbon, meaning carbon molecules with different neutron numbers. The most common is the stable 12C, which make up about 98,9% of the molecules. Out of the non-stable carbon isotopes, a natural production and a long half-life of 14C radiocarbon make the ratio between it and 12C applicable for carbon dating and establishing the amount of fossil fuel burnt, which is the focus here. Carbon is replenished in the tissue as long as the organic object is alive, and that means that the ratio between the two types of carbon isotopes are the essentially the same as the carbon in the air. Once dead and buried with no oxidation, the exchange with the environment stops. The ratio of radiocarbon will then decrease at a known rate – a half-life of 5730 years - and knowing this, the age of the sample can be solved for. In fossil fuel, all the radiocarbon has already decayed and burning it means a shift away from the standard ratio between 12C and 14C in the atmosphere. The shift is measured in per-mil. Measuring the shift from the standard ratio between the two isotopes can in turn be used to calculate how much fossil fuel was burn within the footprint area. However, there are a couple of other factors to consider: The ratio doubled with the atomic bombs, and Chinese bomb tests in 1981, but have since then levelled out. However, there is still an ongoing small decrease with a north-south difference of a few per-mil. These are estimated and added to the standard ratio. Another aspect is the addition of radiocarbon from nuclear power plants and fuel reprocessing stations. The proximity to these will determine their influence, and in this Notebook the influence for aggregated dates with transport models run for specific hours of the day make for an estimated expected shift in permil.
In [1]:
%run ~/STILT_modules_v2.5.ipynb
%run ~/ICOS_atmObs_modules_v2.5.ipynb
#%run ~/ida_master_def_upd.ipynb
%run ~/for_thesis.ipynb
stations = create_STILT_dictionary()
pathFP='/opt/stiltdata/fsicos2/stiltweb/'
from netCDF4 import Dataset
import pandas as pd
import math
import numpy
from mpl_toolkits.basemap import Basemap, cm
import matplotlib
from matplotlib.pyplot import figure
import matplotlib.pyplot as plt
import statistics
In [20]:
list_stations_for_analysis=select_stations()
Once done with station selection, run this cell to finalize selection
It is also possible to write the name of one, or several, more stations to base the analysis on
In [21]:
only_selection_stations=list_stations_for_analysis.value
all_stations=[]
for value in only_selection_stations:
all_stations.append(value)
additional_station=input('Write additional station(s) here, otherwise leave blank. If several, seperate by whitespace: ')
list_additional_station = additional_station.split()
for item in list_additional_station:
all_stations.append(str(item))
print('Your selection is: ', all_stations)
icos_stations = [ist[0:3] for ist in all_stations]
In [4]:
date_range, timeselect, start_date, end_date = date_range_date_hour()
Access and display radiocarbon emission year 2016 based on data from the supplemental material of the paper Zazzeri et al. 2018
"Table S1: emission factor based estimates in TBq/yr calculated by multiplying the standard emission factors in Table 1 by the power production values from the IAEA PRIS Database.
A fraction of 72 % of 14C released as CH4 from the PWRs has been used in the estimation, with all other 14C released in the form of CO2"
Outputs a map with circle sizes corresponding to the amount of emission, and a bar graph showing the names of the emitting sites and their total emissions year 2016.
In [5]:
#update measured emissions for four of the german facilities:
#Gundremmingen: 1073 to 240 GBq
#Isar 2: 87 to 58 Gbq
#Neckarwestheim: 82 to 210 Gbq
#Philippsburg: 89 to 57 GBq
radiocarbon_zazzeri = Dataset('radiocarbon_zazzeri_upd_bfs_2016.nc')
#67 facilities.
lat_lon_site_dic={(47.5416666667, 34.5625): 'Zaporozhe 1-6 (Zaporizhzhia)',(40.9583333333, 0.812499999994): 'Vandellos 2',
(40.7083333333, -2.5625): 'Trillo 1',(44.2916666667, 4.68749999999): 'Tricastin 1-4',( 55.9583333333, -2.43750000001): 'Torness A-B',
(50.5416666667, 5.31249999999): 'Tihange 1-3',(49.2083333333, 14.3125): 'Temelin 1-2',(47.7083333333, 1.56249999999): 'St. Laurent B1-B2',
(45.375, 4.81249999999): 'St. Alban 1-2',(47.7916666667, 31.1875): 'South Ukraine 1-3',
(54.2083333333, 33.1875): 'Smolensk 1',(52.2083333333, 1.56249999999): 'Sizewell B',(54.4583333333, -3.4375): 'Sellafield',(51.2916666667, 25.9375): 'Rovno 1-4 (Rivne)',
(57.2916666667, 12.0625): 'Ringhals 1-4',(49.2916666667, 8.43749999999): 'Philippsburg 2',(49.9583333333, 1.18749999999): 'Penly 1-2',(49.875, 0.687499999994): 'Paluel 1-4',(46.5416666667, 18.8125): 'Paks 1-4',(57.375, 16.6875): 'Oskarshamn 1-3',(61.2083333333, 21.4375): 'Olkiluoto 1-2',(48.5416666667, 3.56249999999): 'Nogent 1-2',(49.0416666667, 9.18749999999): 'Neckarwestheim 2',(46.9583333333, 7.31249999999): 'Muhleberg',(48.2916666667, 18.4375): 'Mochovce',(60.375, 26.3125): 'Loviisa 1-2',(59.875, 29.0625): 'Leningrad 1-4',(47.625, 8.18749999999): 'Leibstadt',(49.7083333333, -1.93750000001): 'La Hague',(45.9583333333, 15.5625): 'Krsko',(43.7083333333, 23.8125): 'Kozloduy 1-6',(67.4583333333, 32.4375): 'Kola 1-4',( 50.2916666667, 26.6875): 'Khmelnitski 1-2',(48.625, 12.3125): 'Isar 2',(55.7083333333, -4.9375): 'Hunterstone B1-B2',(51.2083333333, -3.1875): 'Hinkley Point BA-BB',(54.0416666667, -2.9375): 'Heysham 1A-1B and 2A-2B',(54.625, -1.18750000001): 'Hartlepool A1-A2',(48.5416666667, 10.4375): 'Gundremmingen B-C',(52.0416666667, 9.43749999999): 'Grohnde TOT',(51.0416666667, 2.18749999999): 'Gravelines 1-6',(47.375, 7.93749999999): 'Gosgen',(44.125, 0.812499999994): 'Golfech 1-2',(60.375, 18.1875): 'Forsmark 1-3',( 49.5416666667, -1.93750000001): 'Flamanville 1-2',(47.875, 7.56249999999): 'Fessenheim 1-2',(52.4583333333, 7.31249999999): 'Emsland',(50.875, 0.937499999994): 'Dungeness B1-B2',(49.125, 16.1875): 'Dukovany 1-4',(51.2916666667, 4.31249999999): 'Doel 1-4',(47.7083333333, 2.56249999999): 'Dampierre 1-4',(44.625, 4.81249999999): 'Cruas 1-4',(39.2083333333, -1.06250000001): 'Cofrentes',(46.4583333333, 0.687499999994): 'Civaux',(50.125, 4.81249999999): 'Chooz B1-B2',(47.2083333333, 0.187499999994): 'Chinon B1-B4',(44.2916666667, 28.0625): 'Cernavoda',(49.375, 6.18749999999): 'Cattenom 1-4',( 45.7916666667, 5.31249999999): 'Bugey 2-5',(53.875, 9.31249999999): 'Brokdorf',(51.4583333333, 3.68749999999): 'Borssele',(48.4583333333, 17.6875): 'Bohunice B 1-4', (45.2916666667, -0.687500000006): 'Blayais 1-4',(47.5416666667, 8.18749999999): 'Beznau 1-2',(47.5416666667, 2.81249999999): 'Belleville 1-2',(41.2083333333, 0.562499999994): 'Asco 1-2',( 39.7916666667, -5.6875): 'Almaraz 1-2'}
lon=radiocarbon_zazzeri.variables['lon'][:]
lat=radiocarbon_zazzeri.variables['lat'][:]
fp_radiocarbon = radiocarbon_zazzeri.variables['2016_gigaC'][:,:]
count=0
list_row_col=[]
associated_values=[]
#get the row/col of the sites (only ones that has emission)!
for row in range(fp_radiocarbon.shape[0]):
for col in range(fp_radiocarbon.shape[1]):
if fp_radiocarbon[row][col]>0:
#location of all sources of radiocarbon. for maps.
list_row_col.append([row,col])
#the emission values from the sites.
associated_values.append(fp_radiocarbon[row][col])
#access the latitude and longitude of the cells that have radiocarbon emission in accordace with the above loop.
#it will be one cell per station, hence appropriate to display points for nuclear power plants/reporcessing stations
#using these lat and long values.
latitude_list=[]
longitude_list=[]
for item in list_row_col:
latitude_list.append(lat[item[0]])
longitude_list.append(lon[item[1]])
print("Number of locations: ", len(latitude_list))
#create the map:
figure(num=None, figsize=(18, 16), dpi=80, facecolor='w', edgecolor='k')
map = Basemap(projection='merc', lat_0 = 57, lon_0 = -135,resolution="i",
llcrnrlon=-15, llcrnrlat=33,
urcrnrlon=35, urcrnrlat=73)
map.drawmapboundary(fill_color='lightblue')
map.fillcontinents(color='white')
map.drawcountries()
#sort the data from the cells into differen categories based on emission value (bigger value = show up as bigger circle on map)
less_than_10_latitude_list=[]
less_than_10_longitude_list=[]
less_than_10_associated_values_list=[]
less_than_200_latitude_list=[]
less_than_200_longitude_list=[]
less_than_200_associated_values_list=[]
less_than_400_latitude_list=[]
less_than_400_longitude_list=[]
less_than_400_associated_values_list=[]
less_than_600_latitude_list=[]
less_than_600_longitude_list=[]
less_than_600_associated_values_list=[]
less_than_1000_latitude_list=[]
less_than_1000_longitude_list=[]
less_than_1000_associated_values_list=[]
less_than_2000_latitude_list=[]
less_than_2000_longitude_list=[]
less_than_2000_associated_values_list=[]
more_than_2000_latitude_list=[]
more_than_2000_longitude_list=[]
more_than_2000_associated_values_list=[]
more_than_4000_latitude_list=[]
more_than_4000_longitude_list=[]
more_than_4000_associated_values_list=[]
all_values_list=[]
#list for labels, sting of site name and emission in string.
all_values_list_sites=[]
#maybe add a bar graph too! Put the dictionary up here (lat/long associated with the different stations)
index=0
#the associated values are the values associated with different cells
for value in associated_values:
if value>0:
all_values_list.append(value)
all_values_list_sites.append(str(lat_lon_site_dic[latitude_list[index],longitude_list[index]]) + ' (' + str("%.0f" %value) + ')')
if value<50:
less_than_10_latitude_list.append(latitude_list[index])
less_than_10_longitude_list.append(longitude_list[index])
less_than_10_associated_values_list.append(value)
if value>=50 and value<100:
less_than_200_latitude_list.append(latitude_list[index])
less_than_200_longitude_list.append(longitude_list[index])
less_than_200_associated_values_list.append(value)
if value>=100 and value<200:
less_than_400_latitude_list.append(latitude_list[index])
less_than_400_longitude_list.append(longitude_list[index])
less_than_400_associated_values_list.append(value)
if value>=200 and value<400:
less_than_600_latitude_list.append(latitude_list[index])
less_than_600_longitude_list.append(longitude_list[index])
less_than_600_associated_values_list.append(value)
if value>=400 and value<1000:
less_than_1000_latitude_list.append(latitude_list[index])
less_than_1000_longitude_list.append(longitude_list[index])
less_than_1000_associated_values_list.append(value)
if value>=1000 and value<2000:
less_than_2000_latitude_list.append(latitude_list[index])
less_than_2000_longitude_list.append(longitude_list[index])
less_than_2000_associated_values_list.append(value)
if value>=2000 and value<4000:
more_than_2000_latitude_list.append(latitude_list[index])
more_than_2000_longitude_list.append(longitude_list[index])
more_than_2000_associated_values_list.append(value)
if value>=4000:
more_than_4000_latitude_list.append(latitude_list[index])
more_than_4000_longitude_list.append(longitude_list[index])
more_than_4000_associated_values_list.append(value)
index=index+1
matplotlib.rcParams.update({'font.size': 17})
#different styling dependent on which list the the sites "ended up in"
x_1,y_1 = map(less_than_10_longitude_list, less_than_10_latitude_list)
p1= map.plot(x_1, y_1, 'ro', alpha=0.7, markersize=2)
x_2,y_2 = map(less_than_200_longitude_list, less_than_200_latitude_list)
p2= map.plot(x_2, y_2, 'ro', alpha=0.7, markersize=4)
x_3,y_3 = map(less_than_400_longitude_list, less_than_400_latitude_list)
p3= map.plot(x_3, y_3, "ro", alpha=0.7, markersize=6)
x_4,y_4 = map(less_than_600_longitude_list, less_than_600_latitude_list)
p4= map.plot(x_4, y_4, "ro", alpha=0.7, markersize=8)
x_5,y_5 = map(less_than_1000_longitude_list, less_than_1000_latitude_list)
p5= map.plot(x_5, y_5, "ro", alpha=0.7, markersize=10)
x_6,y_6 = map(less_than_2000_longitude_list, less_than_2000_latitude_list)
p6= map.plot(x_6, y_6, "ro", alpha=0.7, markersize=14)
x_7,y_7 = map(more_than_2000_longitude_list, more_than_2000_latitude_list)
p7= map.plot(x_7, y_7, "ro",alpha=0.7, markersize=18)
x_8,y_8 = map(more_than_4000_longitude_list, more_than_4000_latitude_list)
p8= map.plot(x_8, y_8, "ro",alpha=0.7, markersize=24)
plt.title('Radiocarbon ($^1$$^4$CO$_2$) from nuclear power plants' + '\n' + 'and fuel reprocessing stations, year 2016 (Zazzeri)')
plt.legend((p1[0], p2[0],p3[0], p4[0],p5[0],p6[0],p7[0], p8[0]), ('0 - 50', '50 - 100', '100 - 200', '200 - 400', '400 - 1000', '1000 - 2000', '2000 - 4000' , '19100'), title="Discharge, GBq/Year")
plt.show()
all_values_list_sorted=[x for _,x in sorted(zip(all_values_list,all_values_list))]
all_values_list_sites_sorted=[x for _,x in sorted(zip(all_values_list,all_values_list_sites))]
matplotlib.rcParams.update({'font.size': 17})
x_pos = numpy.arange(len(all_values_list_sorted))
fig, ax = plt.subplots(figsize=(25,10))
ax.bar(x_pos, all_values_list_sorted,
align='center',
alpha=0.5,
color='green')
ax.set_ylabel('GBq/year')
ax.set_xticks(x_pos)
ax.set_xticklabels(all_values_list_sites_sorted, rotation='vertical')
ax.set_title('Radiocarbon ($^1$$^4$CO$_2$) from nuclear power plants' + '\n' + 'and fuel reprocessing stations at stations 2016 (Zazzeri)')
ax.yaxis.grid(True)
plt.show()
Influence of nuclear power plants and fuel reprocessing stations: Time series by station and monthly averages
The emissions from the different sites in the above map need to be translated into what size shift they cause from the standard
ratio of 14C to 12C maintained naturally in the atmosphere. The size of the shift will only depend on the transport to the stations - the footprints - because the emissions do not currently have a temporal resolution.
It might be possible for some stations to get better data in the future.
Fossil fuel emissions are the other anthropogenic influence on the atmospheric 14C to 12C ratio because these do not contain any 14C. Hence, with a measurement of the ratio and establishing the shift from the normal ratio
the fossil fuel emissions within the footprint can be established. Failing to account for the influence from nuclear power plants and fuel reprocessing stations could mean an underestimate.
The emissions are in GBq, and need to be translated into resulting Bq/m3 (activity concentration) they cause at the place of measurement.
$$A_{conc} [\frac{Bq}{m^3}] = Q_{^{14}CO_2} [\frac{Bq}{m^2 s}] \times F \times \frac{1}{V_m} [\frac{mol}{m^3}]$$
within each STILT cell can be calculated using the following
The activity concentration given the standard ratio
The process of translating the emissions in radiocarbon (14CO2) from the facilities listed above consist of first translating the point source emission in
GBq/year into activity concentration in each of the STILT cells given a specific date and time with associated footprint (equation from Ute Karstens):
$$A_{conc} [\frac{Bq}{m^3}] = Q_{^{14}CO_2} [\frac{Bq}{m^2 s}] \times F \times \frac{1}{V_m} [\frac{mol}{m^3}]$$
with molar volume of air (ideal gas) at standard atmospheric pressure and temperature (STP) $V_m = 0.022414 \frac{m^3}{mol}$
The following steps regards translating the emissions into activity concentration in the STILT-cell: $$Q_{^{14}CO_2} [\frac{Bq}{m^2 s}]$$ in the above equation.
GBq --> Bq: multiply by 10^9
per year--> per s: divide by (365 * 24 * 60 * 60)
per point source-->per area of the footprint cell size of the point source: divide by m2 of cell
Once these cell values are known, the activity concentration can be calculated given the transport to the station represented by the footprint for a specific date/time.
The grid with the calculated values is multiplied by the footprint(s) - aggregated and averaged or not - and finally divided by the molar volume of air (0.022414).
Once the activity concentrations in each cell are known, it is possible to calculate $\Delta^{14}C$ using the following equation:
$$\Delta^{14}C = [\frac{Q F}{\chi_{CO_2} M_CA_{abs}}] \times 1000 ‰$$
where
$\chi_{CO_2} [\frac{{\mu}mol_{tracer}}{mol_{air}}]$ $CO_2$ molar mixing ratio. Modelled concentration used generally.
$ M_C [\frac{g}{mol}]$ molar mass of $C$ (12 g)
$V_m [\frac{m^3}{mol}]$ molar Volume of air at standard pressure and temperature STP (0.02447)
$A_{abs} [\frac{Bq}{gC}]$ activity of the standard $A_{abs} = 0.226 [\frac{Bq}{gC}].$ Adjusted to exclude *0.968 in the above equation = 0.238
The questions that are prompted the user when she runs the cell read:
Contribution for date/time to count as a day with influence from nuclear facilities (statistics in graph)
--> the result of this show up in the title of the time series: a percent of how many days have each particular station are influenced by this type of emission.
For some stations a small number of really high influence from these facilities might shift the average, whereas other see this type of influence most date/times.
Average contribution by facility for inclusions (output figure next cell)
--> the time series that are output from running this cell only show the combined influence by all the facilities listed in the bar graph above. However, it is possible to
look at the influence from the separate facilities. The threshold defined here sets a minimum limit to what the average contribution from a facility, during the defined date range,
needs to be to be included in the time series of the next cell.
In [25]:
#threshold for counting a time-step (with assoicated footprint): if the resulting shift is larger than the threshold, added to a count.
#only displayed in the resulting map per station
threshold_days_of_seen=input("Threshold contribution for date/time to count as a day with 'influence from nuclear facilities' (statistics in graph): ")
#in next cell, contribution from the individual facilities. Here threshold defined for being included in the time series.
threshold=input("Average contribution by facility for inclusions (output figure next cell): ")
###################################################################################################################
#Processing to Bq/m2/s in each STILT cell rather than a total GBq of emissions from facilities within the cell.
#GBq--> Bq (total for year 2016)
fp_radiocarbon_bq_yr= fp_radiocarbon*1000000000
#Emission year --> emission per second(365*24*60*60=31536000)
fp_radiocarbon_bq_s=fp_radiocarbon_bq_yr/31536000
#Emission distributed evenly within the STILT cell it is emitted:
f_gridarea = cdf.Dataset('gridareaSTILT.nc')
gridarea = f_gridarea.variables['cell_area'][:]
fp_radiocarbon_bq_s_m2= fp_radiocarbon_bq_s/gridarea
######################################################################################################################
#one list per station with the shift in radiocarbon in permil due to emission from nuclear power plants and fuel reprocessing stations.
list_of_list_delta_radiocarbon=[]
#put all the dates where shift in radiocarbon due to NPPs and fuel reprocessing stations can be calculated
list_of_lists_dates=[]
#same, but for stations (not added to list if no can be calulcated (for instance if STILt has not been run))
list_stations_valid=[]
#for the maps in two cells from now, the values associated with emissions from the different facilites - just need to be over 0:
list_of_lists_row_and_col_associated_val=[]
#the location of the facilities that have emission contribution to shift in radiocarbon at given station for given date range:
list_of_lists_latitude=[]
list_of_lists_longitude=[]
#the total shift per station, given the defined date range
station_average_delta_radiocarbon=[]
station_standard_deviation_radiocarbon=[]
#for time-series for each of the facilties, list of lists of lists because of list per station, then one list for each of the dates,
#and then one list with the values of all the different facilities (given certain date/time)
#associated dates are in the list_of_lists_dates):
list_of_lists_of_lists_row_and_col_associated_val_time_series=[]
list_of_lists_of_lists_row_and_col_time_series=[]
#get the row/col of the facilities that has an average over the threshold
list_of_lists_row_and_col_over_threshold=[]
#Variables used to calculate the expexted shift in radiocarbon
#the volume of air at the standard temperature and pressure
#not used - taken out..
Vm = 0.02447
#what should it be? 0.226? No --> the *0.968 adjusted for (used to noramlize to 13C).
#activity of the standard is 0.226(maintained naturally. How much shift away from this?)
Aabs = 0.238
#molar mass C in gram
Mc = 12
for station in all_stations:
#modelled concentration to calculate the shift in radiocarbon (the natural acitivity dependent on the total C in the air)
modelled_concentration = read_stilt_timeseries(station, date_range)
#get all delta radiocarbon values in a list, one for each station.
#these are the values that goes into each station's time series
list_delta_radiocarbon=[]
#each of the calculated grids are added to this list. Then summed (each cell), added together, and divided by the number of valid dates
#for the station to create an "average grid" --> can get the influence from each facility
list_delta_radiocarbon_contribution_grid=[]
#get a list with the "valid" dates, use in time series.
list_dates=[]
#locations of the facilities that contribute more than the threshold on average
list_row_and_col=[]
#the values of these cells:
list_row_and_col_associated_val=[]
#latitude and longitude can be established based on the row/col values
latitude_list=[]
longitude_list=[]
#for the time series of individual facilities
list_of_lists_row_and_col_associated_val_time_series=[]
list_of_lists_row_and_col_time_series=[]
#list of row/col from the facilities from which we want the values in a time series.
list_row_and_col_over_threshold=[]
#count_days_seen is added to everytime the total contribution from all facilites add up to more than the defined threshold
count_days_seen=0
#want to calculate averages for each of the month. Only add to these if "in the correct month" when looping over the dates
#in the date range. Could do the same for each of the hours.
jan_count=0
jan_tot=0
feb_count=0
feb_tot=0
mar_count=0
mar_tot=0
apr_count=0
apr_tot=0
may_count=0
may_tot=0
jun_count=0
jun_tot=0
jul_count=0
jul_tot=0
aug_count=0
aug_tot=0
sep_count=0
sep_tot=0
oct_count=0
oct_tot=0
nov_count=0
nov_tot=0
dec_count=0
dec_tot=0
#count nan --> number of nan-values by month (will be in labels for the different months)
count_nan_jan=0
count_nan_feb=0
count_nan_mar=0
count_nan_apr=0
count_nan_may=0
count_nan_jun=0
count_nan_jul=0
count_nan_aug=0
count_nan_sep=0
count_nan_oct=0
count_nan_nov=0
count_nan_dec=0
#for calculating the average
number_valid_dates=0
#want to count the nan-values to disply in the time-series
count_nan=0
for date in date_range:
filename=(pathFP+'slots/'+stations[station]['locIdent']+'/'+str(date.year)+'/'+str(date.month).zfill(2)+'/'
+str(date.year)+'x'+str(date.month).zfill(2)+'x'+str(date.day).zfill(2)+'x'+str(date.hour).zfill(2)+'/foot')
#for each date, need lists of contribution from each facility:
#row and col for the faciltiies (as one list item)
list_row_and_col_time_series=[]
#emission from the specific facility
list_row_and_col_associated_val_time_series=[]
if os.path.isfile(filename):
#access the footprint given the date/time. Needed to calculate the influence from the facility.
f_fp = cdf.Dataset(filename)
fp_current=f_fp.variables['foot'][:,:,:]
#get modelled concentrations in correct unit (now in micromole - want it in mole because then use molar mass for C of 12 g)
modelled_concentration_value=modelled_concentration['co2.stilt'][date]/1000000
#keeping it a grid (fp_radiocarbon_bq_s_m2*fp_current) to see the resulting shift in radiocarbon caused by each individual cell
delta_radiocarbon_contribution_grid=(((fp_radiocarbon_bq_s_m2*fp_current) / (modelled_concentration_value* Mc * Aabs)) ) *1000
#get one value summed value of this grid for each date/time. What is used in the time series output from running this cell
delta_radiocarbon_contribution_summed= delta_radiocarbon_contribution_grid[:][:].sum()
#unless it is a nan-value, continue
if math.isnan(delta_radiocarbon_contribution_summed)==False:
#need the dates to correctly show the values in the line graph.
list_dates.append(date)
#needed for display and average (could do a len(list_dates) rather)
number_valid_dates = number_valid_dates + 1
#list_row_col was defined in the last cell --> all the STILT cell with radiocarbon emission within them, their row and col values
#need to do it for each date because want it in time series
for row_col_pair in list_row_col:
#same as "list_row_col"
#take away? replace with list_row_col
list_row_and_col_time_series.append(row_col_pair)
#get the associated value each facility
list_row_and_col_associated_val_time_series.append(delta_radiocarbon_contribution_grid[0][row_col_pair[0]][row_col_pair[1]])
#for averages for the different months
if date.month==1:
jan_count=jan_count+1
jan_tot=jan_tot+delta_radiocarbon_contribution_summed
if date.month==2:
feb_count=feb_count+1
feb_tot=feb_tot+delta_radiocarbon_contribution_summed
if date.month==3:
mar_count=mar_count+1
mar_tot=mar_tot+delta_radiocarbon_contribution_summed
if date.month==4:
apr_count=apr_count+1
apr_tot=apr_tot+delta_radiocarbon_contribution_summed
if date.month==5:
may_count=may_count+1
may_tot=may_tot+delta_radiocarbon_contribution_summed
if date.month==6:
jun_count=jun_count+1
jun_tot=jun_tot+delta_radiocarbon_contribution_summed
if date.month==7:
jul_count=jul_count+1
jul_tot=jul_tot+delta_radiocarbon_contribution_summed
if date.month==8:
aug_count=aug_count+1
aug_tot=aug_tot+delta_radiocarbon_contribution_summed
if date.month==9:
sep_count=sep_count+1
sep_tot=sep_tot+delta_radiocarbon_contribution_summed
if date.month==10:
oct_count=oct_count+1
oct_tot=oct_tot+delta_radiocarbon_contribution_summed
if date.month==11:
nov_count=nov_count+1
nov_tot=nov_tot+delta_radiocarbon_contribution_summed
if date.month==12:
dec_count=dec_count+1
dec_tot=dec_tot+delta_radiocarbon_contribution_summed
#only if not nan...
#this is the values to be displayed in the time series output of this cell
list_delta_radiocarbon.append(delta_radiocarbon_contribution_summed)
#into a list, that then can be summed and divided by the number of valid dates to get an AVERAGE GRID
#which in turn can be used to establish from where the contribution came from (include/exclude facilities for time series and map in next cells).
list_delta_radiocarbon_contribution_grid.append(delta_radiocarbon_contribution_grid)
#for the statistics on how many date/times the station "sees" a significant (user-defined) amout of radiocarbon from emissions of these facilities
if delta_radiocarbon_contribution_summed>float(threshold_days_of_seen):
count_days_seen=count_days_seen+1
#if it is a nan-value (ex will happen if STILT has not been run)
else:
count_nan=count_nan+1
#count nan by month, for graph
if date.month==1:
count_nan_jan=count_nan_jan+1
if date.month==2:
count_nan_feb=count_nan_feb+1
if date.month==3:
count_nan_mar=count_nan_mar+1
if date.month==4:
count_nan_apr=count_nan_apr+1
if date.month==5:
count_nan_may=count_nan_may+1
if date.month==6:
count_nan_jun=count_nan_jun+1
if date.month==7:
count_nan_jul=count_nan_jul+1
if date.month==8:
count_nan_aug=count_nan_aug+1
if date.month==9:
count_nan_sep=count_nan_sep+1
if date.month==10:
count_nan_oct=count_nan_oct+1
if date.month==11:
count_nan_nov=count_nan_nov+1
if date.month==12:
count_nan_dec=count_nan_dec+1
#for each date, append the list of the contribution from the different facilities
list_of_lists_row_and_col_associated_val_time_series.append(list_row_and_col_associated_val_time_series)
#could have appended "list_row_and_col" rather. Just to go with associated values. One whole range of the row- and col values with facilities for each date.
list_of_lists_row_and_col_time_series.append(list_row_and_col_time_series)
#once out of the loop going over the dates, if there are any valid dates continue
if number_valid_dates>0:
#list of stations to associate with the delta radiocarbon values. usually same as stations defined.
list_stations_valid.append(station)
#all the radiocarbon contribution grids added. Result in a summed grid.
summed_radiocarbon_contribution_grid= sum(list_delta_radiocarbon_contribution_grid)
#this is the grid that is looped over later to establish from where the contribution for the given station came from.
average_radiocarbon_contribition_grid=summed_radiocarbon_contribution_grid/number_valid_dates
#this is the average delta radiocarbon value for the station, given the specific date/time range.
summed_average_radiocarbon_contribition_grid=average_radiocarbon_contribition_grid[:][:].sum()
#append it to a list to be accssed in the next cell. No need to list of lists since
#only one value per station. access for title using a station_index.
station_average_delta_radiocarbon.append(summed_average_radiocarbon_contribition_grid)
standard_deviation= numpy.nanstd(list_delta_radiocarbon)
#same as for the average values: no need for a list of lists.
station_standard_deviation_radiocarbon.append(standard_deviation)
#days "seen" --> data/times with a shift larger than the user defined threshold. Here %.
percent_days_seen=(count_days_seen/number_valid_dates)*100
#list_row_col is found in the last cell = NEED to run last cell.
#all the cells where there are a nuclear power plant or fuel reprocessing station.
#this used in two cells from now. Want to map the average influence from facilities.
#checking which facilities should be included in the time series and map (average influce grid checked)... threshold
for row_col_pair in list_row_col:
if average_radiocarbon_contribition_grid[0][row_col_pair[0]][row_col_pair[1]]>0:
#append if over 0
list_row_and_col.append(row_col_pair)
list_row_and_col_associated_val.append(average_radiocarbon_contribition_grid[0][row_col_pair[0]][row_col_pair[1]])
#here, also establosh the coordinate_paris where the contribution is more than defined threshold to see if
#show time series of individual facilities in next cell:
if average_radiocarbon_contribition_grid[0][row_col_pair[0]][row_col_pair[1]]>float(threshold):
list_row_and_col_over_threshold.append(row_col_pair)
#to display on map, need latitude and longitude in seperate lists:
for item in list_row_and_col:
latitude_list.append(lat[item[0]])
longitude_list.append(lon[item[1]])
#create the time series (contribution on station from all facilities combined).
matplotlib.rcParams.update({'font.size': 18})
fig = p.figure(figsize=(20,10))
ax = fig.add_subplot(111)
p.plot(list_dates, list_delta_radiocarbon, linestyle='-',color='r', label='$\Delta^{14}C$ [‰] ')
p.axhline(0, color='black', lw=1)
date_index_number = (len(date_range) - 1)
p.title('Shift in radiocarbon due to emissions at nuclear power plants' + '\n' + 'and fuel reprocessing stations at station ' + str(station) + '\n' + str(date_range[0].year) + '-' + str(date_range[0].month) + '-' + str(date_range[0].day)\
+ ' to ' + str(date_range[date_index_number].year) + '-' + str(date_range[date_index_number].month) + '-' + str(date_range[date_index_number].day)\
+ ', Hour(s): ' + timeselect + '\n' + 'Average: ' + str("%.2f" % summed_average_radiocarbon_contribition_grid) + ' $\Delta^{14}C$ [‰]' + ', Standard deviation: ' + str("%.2f" % standard_deviation) + '\n' + 'Count of days with total average over ' + \
str(threshold_days_of_seen) + ' $\Delta^{14}C$ [‰]: ' + str(count_days_seen) + '/' + str(number_valid_dates) + ' (' + str("%.2f" % percent_days_seen) + '%)' + '\n' + 'Count NaN: ' + str(count_nan))
ax.set_xlim(start_date,end_date)
ax.set_ylabel('$\Delta^{14}C$ [‰] ', color='r')
ax.grid(axis='x')
ax.legend(loc='upper right')
ax.tick_params(axis='y', labelcolor='r')
p.show()
#to use later for maps (where did the contribution come from?). Also to put in
#with the fossil fuel/background data.
list_of_list_delta_radiocarbon.append(list_delta_radiocarbon)
list_of_lists_dates.append(list_dates)
print(len(list_dates))
list_of_lists_row_and_col_associated_val.append(list_row_and_col_associated_val)
list_of_lists_latitude.append(latitude_list)
list_of_lists_longitude.append(longitude_list)
#for the time series showing the contribution by date time from each one of the facilities.
#one final list for each station
list_of_lists_of_lists_row_and_col_associated_val_time_series.append(list_of_lists_row_and_col_associated_val_time_series)
list_of_lists_of_lists_row_and_col_time_series.append(list_of_lists_row_and_col_time_series)
#list over defined threshold
list_of_lists_row_and_col_over_threshold.append(list_row_and_col_over_threshold)
by_month_averages=[]
by_month_names_for_label=[]
if jan_tot>0:
jan_average=jan_tot/jan_count
by_month_averages.append(jan_average)
by_month_names_for_label.append('Jan' + ' (' + str(count_nan_jan) + ')')
if feb_tot>0:
feb_average=feb_tot/feb_count
by_month_averages.append(feb_average)
by_month_names_for_label.append('Feb'+ ' (' + str(count_nan_feb) + ')')
if mar_tot>0:
mar_average=mar_tot/mar_count
by_month_averages.append(mar_average)
by_month_names_for_label.append('Mar'+ ' (' + str(count_nan_mar) + ')')
if apr_tot>0:
apr_average=apr_tot/apr_count
by_month_averages.append(apr_average)
by_month_names_for_label.append('Apr'+ ' (' + str(count_nan_apr) + ')')
if may_tot>0:
may_average=may_tot/may_count
by_month_averages.append(may_average)
by_month_names_for_label.append('May'+ ' (' + str(count_nan_may) + ')')
if jun_tot>0:
jun_average=jun_tot/jun_count
by_month_averages.append(jun_average)
by_month_names_for_label.append('Jun'+ ' (' + str(count_nan_jun) + ')')
if jul_tot>0:
jul_average=jul_tot/jul_count
by_month_averages.append(jul_average)
by_month_names_for_label.append('Jul'+ ' (' + str(count_nan_jul) + ')')
if aug_tot>0:
aug_average=aug_tot/aug_count
by_month_averages.append(aug_average)
by_month_names_for_label.append('Aug' + ' (' + str(count_nan_aug) + ')')
if sep_tot>0:
sep_average=sep_tot/sep_count
by_month_averages.append(sep_average)
by_month_names_for_label.append('Sep'+ ' (' + str(count_nan_sep) + ')')
if oct_tot>0:
oct_average=oct_tot/oct_count
by_month_averages.append(oct_average)
by_month_names_for_label.append('Oct'+ ' (' + str(count_nan_oct) + ')')
if nov_tot>0:
nov_average=nov_tot/nov_count
by_month_averages.append(nov_average)
by_month_names_for_label.append('Nov'+ ' (' + str(count_nan_nov) + ')')
if dec_tot>0:
dec_average=dec_tot/dec_count
by_month_averages.append(dec_average)
by_month_names_for_label.append('Dec' + ' (' + str(count_nan_dec) + ')')
N = len(by_month_averages)
ind = np.arange(N)
width = 0.35
from matplotlib.pyplot import figure
figure(num=None, figsize=(11, 7), dpi=80, facecolor='w', edgecolor='k')
ax = plt.subplot(111)
ax.yaxis.grid(True)
plt.bar(ind, by_month_averages, width, color='blue')
matplotlib.rcParams.update({'font.size': 12})
plt.ylabel('$\Delta^{14}C$ [‰]')
plt.xlabel('Month (count NaN)')
#added figsize
#for title --> use variables of last cell as well (ex number of no measurements. nan here can also just mean no cells for it)
plt.title(station +' average influence by month' + '\n' + str(date_range[0].year) + '-' + str(date_range[0].month) + '-' + str(date_range[0].day)\
+ ' to ' + str(date_range[date_index_number].year) + '-' + str(date_range[date_index_number].month) + '-' + str(date_range[date_index_number].day) + ', Hour(s): ' + timeselect)
plt.xticks(ind, by_month_names_for_label, rotation='vertical')
plt.axhline(0, color='black', lw=2)
plt.show()
Influence of nuclear power plants and fuel reprocessing stations: Time series by station showing individual facilities
The values for this cell come from running the previous cell. That is also where the threshold for inclusion in the time series is defined. Note, if you want to change the threshold you need to run the last cell one more time.
In [27]:
#dictionary to use the lat and lon values of each of the facilities to access their names.
lat_lon_site_dic={(47.5416666667, 34.5625): 'Zaporozhe 1-6 (Zaporizhzhia)',(40.9583333333, 0.812499999994): 'Vandellos 2',(40.7083333333, -2.5625): 'Trillo 1',(44.2916666667, 4.68749999999): 'Tricastin 1-4',( 55.9583333333, -2.43750000001): 'Torness A-B',(50.5416666667, 5.31249999999): 'Tihange 1-3',(49.2083333333, 14.3125): 'Temelin 1-2',(47.7083333333, 1.56249999999): 'St. Laurent B1-B2',(45.375, 4.81249999999): 'St. Alban 1-2',(47.7916666667, 31.1875): 'South Ukraine 1-3',(54.2083333333, 33.1875): 'Smolensk 1',(52.2083333333, 1.56249999999): 'Sizewell B',(54.4583333333, -3.4375): 'Sellafield',(51.2916666667, 25.9375): 'Rovno 1-4 (Rivne)',(57.2916666667, 12.0625): 'Ringhals 1-4',(49.2916666667, 8.43749999999): 'Philippsburg 2',(49.9583333333, 1.18749999999): 'Penly 1-2',(49.875, 0.687499999994): 'Paluel 1-4',(46.5416666667, 18.8125): 'Paks 1-4',(57.375, 16.6875): 'Oskarshamn 1-3',(61.2083333333, 21.4375): 'Olkiluoto 1-2',(48.5416666667, 3.56249999999): 'Nogent 1-2',(49.0416666667, 9.18749999999): 'Neckarwestheim 2',(46.9583333333, 7.31249999999): 'Muhleberg',(48.2916666667, 18.4375): 'Mochovce',(60.375, 26.3125): 'Loviisa 1-2',(59.875, 29.0625): 'Leningrad 1-4',(47.625, 8.18749999999): 'Leibstadt',(49.7083333333, -1.93750000001): 'La Hague',(45.9583333333, 15.5625): 'Krsko',(43.7083333333, 23.8125): 'Kozloduy 1-6',(67.4583333333, 32.4375): 'Kola 1-4',( 50.2916666667, 26.6875): 'Khmelnitski 1-2',(48.625, 12.3125): 'Isar 2',(55.7083333333, -4.9375): 'Hunterstone B1-B2',(51.2083333333, -3.1875): 'Hinkley Point BA-BB',(54.0416666667, -2.9375): 'Heysham 1A-1B and 2A-2B',(54.625, -1.18750000001): 'Hartlepool A1-A2',(48.5416666667, 10.4375): 'Gundremmingen B-C',(52.0416666667, 9.43749999999): 'Grohnde TOT',(51.0416666667, 2.18749999999): 'Gravelines 1-6',(47.375, 7.93749999999): 'Gosgen',(44.125, 0.812499999994): 'Golfech 1-2',(60.375, 18.1875): 'Forsmark 1-3',( 49.5416666667, -1.93750000001): 'Flamanville 1-2',(47.875, 7.56249999999): 'Fessenheim 1-2',(52.4583333333, 7.31249999999): 'Emsland',(50.875, 0.937499999994): 'Dungeness B1-B2',(49.125, 16.1875): 'Dukovany 1-4',(51.2916666667, 4.31249999999): 'Doel 1-4',(47.7083333333, 2.56249999999): 'Dampierre 1-4',(44.625, 4.81249999999): 'Cruas 1-4',(39.2083333333, -1.06250000001): 'Cofrentes',(46.4583333333, 0.687499999994): 'Civaux',(50.125, 4.81249999999): 'Chooz B1-B2',(47.2083333333, 0.187499999994): 'Chinon B1-B4',(44.2916666667, 28.0625): 'Cernavoda',(49.375, 6.18749999999): 'Cattenom 1-4',( 45.7916666667, 5.31249999999): 'Bugey 2-5',(53.875, 9.31249999999): 'Brokdorf',(51.4583333333, 3.68749999999): 'Borssele',(48.4583333333, 17.6875): 'Bohunice B 1-4', (45.2916666667, -0.687500000006): 'Blayais 1-4',(47.5416666667, 8.18749999999): 'Beznau 1-2',(47.5416666667, 2.81249999999): 'Belleville 1-2',(41.2083333333, 0.562499999994): 'Asco 1-2',( 39.7916666667, -5.6875): 'Almaraz 1-2'}
station_index=0
for station in all_stations:
#get the facilities that has an average over the user-defined threshold
facility_list=[]
for coordinate_pair in list_of_lists_row_and_col_over_threshold[station_index]:
facility_list.append(lat_lon_site_dic[lat[coordinate_pair[0]],lon[coordinate_pair[1]]])
date_index=0
#for each date/time (footprint), there will be one value per facility. List of list with these:
list_of_lists_values_for_time_series=[]
#for every date/time, a list of contribution from all different facilties.
for date in list_of_lists_dates[station_index]:
#one list for each date/time (footprint)
values_for_times_series_by_date=[]
#row_col --> latitude and longitude associated with the facility locations that have an higher average
#than the user-defined threshold.
for row_col in list_of_lists_row_and_col_over_threshold[station_index]:
if row_col in list_of_lists_of_lists_row_and_col_time_series[station_index][date_index]:
#print(list_of_lists_of_lists_row_and_col_time_series[station_index][date_index])
index_within=list_of_lists_of_lists_row_and_col_time_series[station_index][date_index].index(row_col)
values_for_times_series_by_date.append(list_of_lists_of_lists_row_and_col_associated_val_time_series[station_index][date_index][index_within])
date_index=date_index+1
list_of_lists_values_for_time_series.append(values_for_times_series_by_date)
#only draw a time series with line for each facility if there are any facilities wiht an average contribution higher than the
#user-defined threshold
if len(facility_list)>0:
matplotlib.rcParams.update({'font.size': 17})
fig = p.figure(figsize=(20,10))
ax = fig.add_subplot(111)
facility_index=0
for facility in facility_list:
list_by_facility=[]
for list_item in list_of_lists_values_for_time_series:
if not list_item:
list_by_facility.append(0)
else:
list_by_facility.append(list_item[facility_index])
p.plot(list_of_lists_dates[station_index],
list_by_facility, linestyle='-',
color=numpy.random.rand(3,),
label=(str(facility) + ', average: ' + str("%.4f" %numpy.nanmean(list_by_facility)) + ' $\Delta^{14}C$ [‰]'))
facility_index=facility_index+1
p.axhline(0, color='black', lw=1)
date_index_number = (len(date_range) - 1)
p.title('Shift in radiocarbon due to emissions at nuclear power plants' + '\n' + 'and fuel reprocessing stations at station ' + str(station) + '\n' + str(date_range[0].year) + '-' + str(date_range[0].month) + '-' + str(date_range[0].day)\
+ ' to ' + str(date_range[date_index_number].year) + '-' + str(date_range[date_index_number].month) + '-' + str(date_range[date_index_number].day)\
+ ', Hour(s): ' + timeselect + '\n' + 'Average by facility for inclusion: ' + str(threshold) + ' $\Delta^{14}C$ [‰] ' + '\n' + 'Count NaN: ' + str(count_nan))
ax.set_xlim(start_date,end_date)
ax.set_ylabel('$\Delta^{14}C$ [‰] ')
#ax.grid(axis='y')
ax.legend(loc='upper left')
ax.tick_params(axis='y')
ax.yaxis.grid(True)
p.show()
#if a station has no facilities with an average over the user-defined threshold, give them this message:
else:
print('No facilities at ' + str(station) + ' with an average higher than ' + str(threshold) + ' permil')
#added to every time new station in all_stations looped over.
#needed to move to the next longitude/latitude lists
station_index=station_index+1
Average shift in radiocarbon due to emission at nuclear power plants and fuel reprocessing stations: Maps and bar graphs by station.
Need to run above cells. The values in the bar graphs are the same as the averages by facilities in last the cell.
In [28]:
#the latitude and longitude values of stations - for points on the map.
stations_lat=[]
stations_lon=[]
#access the lat and long for all VALID stations (that have footprints in the defined date range, determined in the last section)
#"original" refers to the height of the station(HTM150 as opposed to shorter HTM)
#want to show the station on the final maps.
for station in list_stations_valid:
stations_lat.append(stations[station]['lat'])
stations_lon.append(stations[station]['lon'])
#to access move in the list
index_list=0
#one list for each station = no need loop over the stations as well (list of stations used to access station "we are on" in the loop)
for list_item in list_of_lists_latitude:
latitude_list=list_of_lists_latitude[index_list]
longitude_list=list_of_lists_longitude[index_list]
#a whole list with values - all over 0.
associated_values=list_of_lists_row_and_col_associated_val[index_list]
#access specific lat/long for station.
station_lat=stations_lat[index_list]
station_lon=stations_lon[index_list]
#here access correct station.
station=list_stations_valid[index_list]
#average radiocarbon contribution from nuclear power plants and fuel reprocessing stations at specific location for soecified date/time range
total_station=station_average_delta_radiocarbon[index_list]
standard_deviation_station= station_standard_deviation_radiocarbon[index_list]
index_list=index_list+1
figure(num=None, figsize=(20, 15), dpi=80, facecolor='w', edgecolor='k')
map = Basemap(projection='merc', lat_0 = 57, lon_0 = -135,resolution="i",
llcrnrlon=-15, llcrnrlat=33,
urcrnrlon=35, urcrnrlat=73)
map.drawmapboundary(fill_color='lightblue')
map.fillcontinents(color='white')
map.drawcountries()
step_1_latitude_list=[]
step_1_longitude_list=[]
step_1_associated_values_list=[]
step_2_latitude_list=[]
step_2_longitude_list=[]
step_2_associated_values_list=[]
step_3_latitude_list=[]
step_3_longitude_list=[]
step_3_associated_values_list=[]
step_4_latitude_list=[]
step_4_longitude_list=[]
step_4_associated_values_list=[]
step_5_latitude_list=[]
step_5_longitude_list=[]
step_5_associated_values_list=[]
step_6_latitude_list=[]
step_6_longitude_list=[]
step_6_associated_values_list=[]
over_0_05_list=[]
over_0_05_list_sites=[]
over_0_05_list_percent=[]
index=0
for value in associated_values:
if value>0.05:
over_0_05_list.append(value)
over_0_05_list_percent.append((value/total_station)*100)
over_0_05_list_sites.append(str(lat_lon_site_dic[latitude_list[index],longitude_list[index]]))
if value>0.05 and value<=0.2:
step_1_latitude_list.append(latitude_list[index])
step_1_longitude_list.append(longitude_list[index])
step_1_associated_values_list.append(value)
if value>0.2 and value<=0.3:
step_2_latitude_list.append(latitude_list[index])
step_2_longitude_list.append(longitude_list[index])
step_2_associated_values_list.append(value)
if value>0.3 and value<=0.4:
step_3_latitude_list.append(latitude_list[index])
step_3_longitude_list.append(longitude_list[index])
step_3_associated_values_list.append(value)
if value>0.4 and value<0.5:
step_4_latitude_list.append(latitude_list[index])
step_4_longitude_list.append(longitude_list[index])
step_4_associated_values_list.append(value)
if value>0.5 and value<=0.6:
step_5_latitude_list.append(latitude_list[index])
step_5_longitude_list.append(longitude_list[index])
step_5_associated_values_list.append(value)
if value>0.6:
step_6_latitude_list.append(latitude_list[index])
step_6_longitude_list.append(longitude_list[index])
step_6_associated_values_list.append(value)
index=index+1
x_1,y_1 = map(step_1_longitude_list, step_1_latitude_list)
p1= map.plot(x_1, y_1, 'ro', alpha=0.7, markersize=3)
x_2,y_2 = map(step_2_longitude_list, step_2_latitude_list)
p2= map.plot(x_2, y_2, 'ro', alpha=0.7, markersize=5)
x_3,y_3 = map(step_3_longitude_list, step_3_latitude_list)
p3= map.plot(x_3, y_3, "ro",alpha=0.7, markersize=7)
x_4,y_4 = map(step_4_longitude_list, step_4_latitude_list)
p4= map.plot(x_4, y_4, "ro",alpha=0.7, markersize=10)
x_5,y_5 = map(step_5_longitude_list, step_5_latitude_list)
p5= map.plot(x_5, y_5, "ro",alpha=0.7, markersize=14)
x_6,y_6 = map(step_6_longitude_list, step_6_latitude_list)
p6= map.plot(x_6, y_6, "ro",alpha=0.7, markersize=18)
x_7,y_7 = map(station_lon, station_lat)
p7= map.plot(x_7, y_7, "yo", markersize=11, alpha=.8)
plt.title('Shift in radiocarbon at ' +str(station) + '\n' + 'due to emissions at nuclear power plants and fuel reprocessing stations' + '\n' + str(date_range[0].year) + '-' + str(date_range[0].month) + '-' + str(date_range[0].day)\
+ ' to ' + str(date_range[date_index_number].year) + '-' + str(date_range[date_index_number].month) + '-' + str(date_range[date_index_number].day)\
+ ', Hour(s): ' + timeselect + '\n' + 'Total average $\Delta^{14}C$ [‰]: ' + str("%.2f" % total_station) + ', Standard deviation: ' + str("%.2f" % standard_deviation_station))
#come back and fix this.
#plt.legend('dots corresponding to amout emitted')
plt.legend((p1[0], p2[0],p3[0], p4[0],p5[0], p6[0], p7[0]), ('0.05 - 0.2', '0.2 - 0.3', '0.3 - 0.4', '0.4 - 0.5', '0.5 - 0.6', '> 0.6', str(station)), title="$\Delta^{14}C$ [‰]")
plt.show()
over_0_05_list_sorted = [x for _,x in sorted(zip(over_0_05_list,over_0_05_list))]
over_0_05_list_sites_sorted = [x for _,x in sorted(zip(over_0_05_list,over_0_05_list_sites))]
over_0_05_list_percent_sorted = [x for _,x in sorted(zip(over_0_05_list,over_0_05_list_percent))]
matplotlib.rcParams.update({'font.size': 17})
x_pos = numpy.arange(len(over_0_05_list))
fig, ax = plt.subplots(figsize=(20,10))
ax.bar(x_pos, over_0_05_list_sorted,
align='center',
alpha=0.5,
color='red')
ax.set_ylabel('$\Delta^{14}C$ [‰]')
ax.set_xticks(x_pos)
ax.set_xticklabels(over_0_05_list_sites_sorted, rotation='vertical')
ax.set_title('Shift in radiocarbon at ' +str(station) + '\n' + 'due to emissions at nuclear power plants and fuel reprocessing stations' + '\n' + str(date_range[0].year) + '-' + str(date_range[0].month) + '-' + str(date_range[0].day)\
+ ' to ' + str(date_range[date_index_number].year) + '-' + str(date_range[date_index_number].month) + '-' + str(date_range[date_index_number].day) + \
', Hour(s): ' + timeselect + '\n' + ' Total average $\Delta^{14}C$ [‰]: ' + str("%.2f" % total_station) + '\n' + '(stations with contribution >0.05 included)')
ax.yaxis.grid(True)
#instantiate a second axes that shares the same x-axis
ax2 = ax.twinx()
ax2.set_ylabel('Percent of total contribution')
ax2.bar(x_pos, over_0_05_list_percent_sorted,color='red',label='Percent')
ax2.tick_params(axis='y')
plt.show()
