calculate hazard curves from GMFs

[Muharem Hrnjadovic]

Implement Hazard curve calculation from GMFs in Event Based calculator

[Damiano Monelli (monelli)]

The extraction of hazard curves from ground motion fields (gmfs) should be an option that a user can select in the [output] section of the job configuration file. That is a user should be able to specify (as in the classical calculation mode):

intensity_measure_types_and_levels = {"PGA": [0.005, 0.007, 0.0098, 0.0137, 0.0192, 0.0269, 0.0376, 0.0527, 0.0738, 0.103, 0.145, 0.203, 0.284, 0.397, 0.556, 0.778, 1.09, 1.52, 2.13], "SA(0.025)": [0.005, 0.007, 0.0098, 0.0137, 0.0192, 0.0269, 0.0376, 0.0527, 0.0738, 0.103, 0.145, 0.203, 0.284, 0.397, 0.556, 0.778, 1.09, 1.52, 2.13]}

if intensity_measure_types_and_levels is empty than this means that hazard curves are not extracted. There must be consistency between the intensity measure types defined in intensity_measure_types_and_levels and in intensity_measure_types (in the [calculation] section). That is if gmfs are computed only for PGA, a user cannot extract the hazard curves for SA.

By default hazard curves should be computed for all the logic tree paths and also from the complete logic tree gmf set (that is the set of gmfs obtained by joining the gmfs from all logic tree paths)

The extraction of hazard curves from a set of gmf is descrived below.

Given a site, an intensity measure type (PGA or SA for instance), a list of intensity measure levels ([IML_1, IML_2, ..., IML_N]) extract all the computed ground motion values ([GMV_1, GMV_2, ..., GMV_M]) and the corresponding duration T (if the ground motion fields correspond to a logic tree path, the duration is 'ses_per_logic_tree_path * investigation_time', if the gmf set corresponds to the complete logic tree, the duration is 'ses_per_logic_tree_path * investigation_time * number_of_logic_tree_samples' if Monte Carlo sampling of logic tree is chosen or 'ses_per_logic_tree_path * investigation_time * number_of_logic_tree_paths' if Enumeration is chosen)

For each IML in [IML_1, IML_2, ..., IML_N] compute the number of values (N) in [GMV_1, GMV_2, ..., GMV_M] that are greater or equal than IML. The probability of exceeding IML can be then computed as:
P(>=IML) = 1 - exp(- (investigation_time / duration) * N )

[Damiano Monelli (monelli)]

Test data for hazard curve calculation with ground motion fields.

The files gmf_values_site1.dat and gmf_values_site2.dat contains ground motion field values for 2 sites both associated to a duration T = 1000 year.

For intensity measure levels [0.01, 0.1, 0.2] the probabilities of exceedance for site 1 and investigation time equal to 1 year are (rounding to 5 digits):
[0.63578, 0.39347, 0.07965]

and for site 2:
[0.63578, 0.28609, 0.02664]

[Lars Butler (lars-butler)]

First pull request:
https://github.com/gem/oq-engine/pull/877

[Lars Butler (lars-butler)]

Had a discussion with al-maisan and monelli about the size of the data sets we'll have to work with. Here's what we discussed:

For a realistic event-based calculation, we can expect
  * A few 100,000s of points
  * For each point, a few 100,000s to a few 1,000,000s of GMVs

Thus, fetching this massive amount of data from the DB, number crunching in Python land, and writing back to the DB perhaps is not the ideal solution.

al-maisan and I discussed this and I'll experiment a bit with doing this number crunching in the database (as a plpython function, for example).

[Lars Butler (lars-butler)]

Slightly tangential, but related: https://github.com/gem/oq-engine/pull/886

[Lars Butler (lars-butler)]

Last patch: https://github.com/gem/oq-engine/pull/913