Skip to content

Basics: building blocks and recipes

This page collects small, self-contained recipes for the individual components of the package. Every code block on this page is executed by the test-suite, so they can be copied as is. The blocks of a section build on the previous ones, as in a notebook.

For complete scenarios, see the short term and long term pages.

Installation

The package targets Python 3.9. From a clone of the repository:

git clone https://github.com/duxfrederic/ensi_g14_calculator.git
cd ensi_g14_calculator
python -m pip install .

or directly from GitHub:

python -m pip install git+https://github.com/duxfrederic/ensi_g14_calculator.git

Check that it works:

import ensi_g14_calculator
from ensi_g14_calculator.short_term_scenarios.runner import ScenarioEvaluator
print("import ok")

Units

Every physical quantity passed to the package carries an astropy unit; results come back with units as well. The sievert is not available in the astropy version supported (Python 3.9), so the package defines it: import it with from ensi_g14_calculator.misc import Sv.

import numpy as np
from astropy.units import Bq, m, s, mm, year
from ensi_g14_calculator.misc import Sv

activity = 1e9 * Bq
distance = 500 * m
print(distance.to('km'))             # 0.5 km, any length unit is converted internally

dose = 3.2e-6 * Sv
print(dose.value)                    # 3.2e-06, the bare float
print(dose.to_value(Sv) * 1e3)       # 0.0032, in mSv (the prefixed units are not registered by name, scale by hand)

A short term dose at a single point

The ScenarioEvaluator accepts scalars, lines or grids for x and y (positions in metres, east and north of the chimney base).

from ensi_g14_calculator.orography import Orography
from ensi_g14_calculator.averaging.weather_averager import SourceParams
from ensi_g14_calculator.short_term_scenarios.runner import ScenarioEvaluator

source = SourceParams(
    H_a=75 * m,       # height of the chimney (total of the stack)
    H_G=15 * m,       # height of the tallest adjacent buildings or structures
    D_i=1.2 * m,      # internal diameter of the chimney
    D_e=1.2 * m,      # external diameter (if unknown, take D_e = D_i as per the directive)
    W_0=10.0 * m / s, # exhaust speed of the chimney
    # z_m=10 * m,     # optional: wind measurement height, default 10 m
)
evaluator = ScenarioEvaluator(source, Orography(flat=True))

nuclides = {'Cs-137': 1e10 * Bq}
results = evaluator.evaluate_short_term_dose(
    nuclides_bq=nuclides,
    x=800.0 * m,          # 800 m east of the chimney ...
    y=0.0 * m,
    wind_direction=-90,   # ... and the wind blows from the west, i.e. towards the point.
)

dose_first_8h = results.sum_doses(weather='D', interval='0-8h', age_group='adult')
print(dose_first_8h)

wind_direction is the direction the wind comes from, in degrees (north = 0, east = 90); the default -90 is a wind from the west, so a point to the east of the chimney is downwind.

Where does the dose come from?

sum_doses takes optional filters, so you can split the dose by pathway, nuclide, interval or age group. The pathways are returned by evaluator.short_term_dose; the safest way to list what a result holds is to look at its keys:

keys = list(results._doses)           # DoseKey(weather, nuclide, interval, pathway, age_group, model, is_reference)
pathways = sorted({k.pathway for k in keys})
intervals = sorted({k.interval for k in keys})
print(pathways)
print(intervals)
print(evaluator.age_groups)

for pathway in pathways:
    d = results.sum_doses(weather='D', interval='0-8h', age_group='adult', pathway=pathway)
    if d is not None:
        print(f"{pathway:>16s}: {float(d.to_value(Sv)):.3e} Sv")

sum_doses returns None when no record matches the filters (e.g. a pathway that does not exist for the requested interval).

Comparing the weather categories and the age groups

for weather in 'ABCDEF':
    adult = results.sum_doses(weather=weather, interval='0-8h', age_group='adult')
    toddler = results.sum_doses(weather=weather, interval='0-8h', age_group='toddler_1_year')
    print(weather, f"adult: {float(adult.to_value(Sv)):.2e} Sv", f"toddler: {float(toddler.to_value(Sv)):.2e} Sv")

The age groups are 'adult', 'child_10_year' and 'toddler_1_year'.

A dose profile along a line

Passing arrays gives a profile; here the first-8h dose along the downwind axis:

import matplotlib
matplotlib.use("Agg")   # no display needed
import matplotlib.pyplot as plt

x = np.linspace(200, 5000, 100) * m
y = np.zeros(100) * m
profile = evaluator.evaluate_short_term_dose(nuclides_bq=nuclides, x=x, y=y, wind_direction=-90)

fig, ax = plt.subplots()
for weather in 'ABCDEF':
    ax.loglog(x.to_value(m), profile.sum_doses(weather=weather, interval='0-8h', age_group='adult').to_value(Sv),
              label=weather)
ax.set_xlabel('distance downwind [m]')
ax.set_ylabel('dose, first 8 h [Sv]')
ax.legend()
fig.savefig('profile.png')

Changing the exposure time and the release date

The exposure time T_exp must be exactly 1 year (default) or 50 years, as in the directive. The release date is a 'mm-dd' string (default '08-01') and determines which crops are exposed.

results_50y = evaluator.evaluate_short_term_dose(
    nuclides_bq=nuclides, x=800.0 * m, y=0.0 * m,
    T_exp=50 * year,
    release_date='06-15',
)
print(results_50y.sum_doses(weather='D', interval='medium_term', age_group='adult'))

Using the wind speed override

The directive prescribes a wind of 1 m/s for short term scenarios. Should you need to deviate, set it on the evaluator (you are then leaving the scope of the directive):

custom = ScenarioEvaluator(source, Orography(flat=True))
custom.wind_speed = 3.0 * m / s
res = custom.evaluate_short_term_dose(nuclides_bq=nuclides, x=800.0 * m, y=0.0 * m)
print(res.sum_doses(weather='D', interval='0-8h', age_group='adult'))

Geometry: coordinates and elevation

Orography converts between WGS84 coordinates and the local metric frame centred on the chimney, and interpolates elevation from an ASC file (see the long term page for how to get one). The block below uses a small elevation map around Bern, shipped in the test data; adapt the path and coordinates to your own file.

chimney_lat, chimney_lon = 46.93026329420960, 7.420436181732923
oro = Orography(asc_path='data/orography.asc', latitude0=chimney_lat, longitude0=chimney_lon)

group_lat, group_lon = 46.92626527964869, 7.41639267991673
gx, gy = oro.latlon_to_xy(group_lat, group_lon)
print(f"group at x={gx:.0f}, y={gy:.0f}")
print("elevation relative to chimney base:", oro.elevation_xy(gx, gy))
print("inside raster:", oro.contains_latlon(group_lat, group_lon))

Evaluate a short term scenario for that real location. Because the orography is used, the terrain influences the dispersion:

terrain_evaluator = ScenarioEvaluator(source, oro)
at_group = terrain_evaluator.evaluate_short_term_dose(
    nuclides_bq=nuclides, x=gx, y=gy, wind_direction=45.0
)
print(at_group.sum_doses(weather='D', interval='0-8h', age_group='adult'))

Elevation profile between two points:

distance_along, elevation = oro.profile_along_line_xy(0 * m, 0 * m, gx, gy, n=50)
print(distance_along[-1], elevation.max())

The dispersion factor on its own

If you only need the atmospheric dispersion \(\chi\) (\(\mathrm{s/m^3}\)), use the WeatherAverager with a weather statistic given as a table. Here a single bin (category D, 3 m/s, wind from the west):

import pandas as pd
from ensi_g14_calculator.averaging.weather_averager import WeatherAverager, WeatherStats

stats = WeatherStats.from_dataframe(
    pd.DataFrame({
        'weather_category': ['D'],
        'wind_speed': [3.0],
        'wind_direction': [270.0],
        'probability': [1.0],
    }),
    randomize_wind_direction_within_bins=False,
)
averager = WeatherAverager(orography=Orography(flat=True), source_params=source)

xg, yg = np.meshgrid(np.linspace(-1000, 1000, 5), np.linspace(-1000, 1000, 5))
chi_k, chi_ks = averager.average_dispersion(xg * m, yg * m, stats)
print(chi_k.unit, chi_k.shape)
print(chi_k.max())

chi_k is the dispersion factor and chi_ks the pseudo-dispersion used for the submersion pathway. With several weather bins the result is the probability-weighted average, which is what the long term scenarios use.

Reading the results of a long term run

After a LongTermEvaluator.run (see the long term page), everything is in an sqlite3 file. For an interactive breakdown of the dose (by chimney, nuclide, pathway, ...), use the dashboards in ensi_g14_dashboard, see its readme.

To script your own analyses instead, the file can be queried from Python. Here we run a small case with the test data, then query it:

import pandas as pd
from astropy.units import m
from ensi_g14_calculator.long_term_scenarios import LongTermEvaluator

weather = pd.read_csv('data/weather_statistics.csv')
LE = LongTermEvaluator('results.sqlite3', 'data/orography.asc')
LE.run('data/chimneys.csv', 'data/exposed_groups.csv',
       weather.drop('rainfall_intensity', axis=1), weather[weather['rainfall_intensity'] > 0.0],
       'data/emission.csv', summer_total_precipitation=0.3 * m)

from ensi_g14_calculator.dose.dose_database import DoseDatabase

with DoseDatabase('results.sqlite3') as db:
    print(db.get_total_dose('mcdo', 'adult'))            # all chimneys / nuclides / pathways
    print(db.get_dose_by_pathway('mcdo', 'adult'))       # ... split by pathway
    print(db.get_dose_by_nuclide('mcdo', 'adult'))       # ... by nuclide
    print(db.get_dose_by_chimney('mcdo', 'adult'))       # ... by chimney
    print(db.get_dose_breakdown('mcdo', 'adult', nuclide='H-3_HTO'))   # one row per chimney/nuclide/pathway

All the query methods return pandas dataframes, with doses in Sv. Note that the pathway names of the long term database (e.g. ingestion_milk_dairy) are not identical to the short term ones (ingestion_dairy). The same table can be read directly with sqlite3 or pandas:

import sqlite3
with sqlite3.connect('results.sqlite3') as conn:
    doses = pd.read_sql('SELECT * FROM doses', conn)
print(doses.columns.tolist())
print(doses.groupby(['group_id', 'age_group'])['dose_Sv'].sum())