vsm/plot_fmd.py at main · guyomd/vsm · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
import os
import glob
import numpy as np
from shapely import Point, Polygon, MultiPolygon, unary_union
from argparse import ArgumentParser

# Internal dependencies:
from lib.plotutils import fmd_histogram
from lib.ioutils import ParameterSet, load_points, load_bins, minmax_in_polygon_file, load_polygons
from lib.geoutils import convert_to_EPSG
from compute_ab_values import TruncatedGRestimator


def _set_value_if_not_nan(value, additional_condition=True):
    if np.isnan(value):
        return None
    elif additional_condition:
        return value


if __name__ == "__main__":

    # Read input arguments:
    parser = ArgumentParser(description="Plot Frequency-Magnitude Distributions (FMD) for a set of pixels in a zoneless model")
    parser.add_argument("configfile",
                        help="Configuration file")

    parser.add_argument("index",
                        help="Cell indices (starting index at 1)",
                        nargs='+',
                        type=int)

    parser.add_argument("-g", "--showgrid",
                        help="Show grid",
                        action="store_true")

    parser.add_argument('-o', '--output-directory',
                        help="Output directory for figures",
                        default=None)

    parser.add_argument('--no-model',
                        help="Do not overlay FMD model adjustment",
                        action='store_true')

    parser.add_argument("-u", "--use-bootstrap-results",
                        help="Read a- and b-values from the mixture distribution or from bootstrapped results",
                        nargs=1,
                        default=None,
                        choices=['mixture', 'bootstrap'])

    args = parser.parse_args()

    # Load parameters:
    prms = ParameterSet()
    prms.load_settings(args.configfile)

    # Check output directory existence:
    if args.output_directory is None:
        args.output_directory = prms.figures_dir
    if not os.path.isdir(args.output_directory):
        os.mkdir(args.output_directory)

    # Load geographical boundaries:
    if prms.bounds_file is not None:
        # Regular grid case:
        envelope = load_points(prms.bounds_file)  # returns a shapely.Polygon instance
    else:
        # Polygonal cells case:
        multipol, _ = load_polygons(prms.mesh_file)
        envelope = unary_union(multipol)
        if isinstance(envelope, MultiPolygon):
            raise Warning(f'Disjoint polygonal cells in file "{prms.mesh_file}". Please fix this.')


    # Load seismicity rates and truncated G-R parameters:
    inputfile = os.path.join(prms.output_dir, 'gridded_densities.txt')
    print(f'>> INPUT: Load densities from file {inputfile}')
    estim = TruncatedGRestimator()
    # --> Rescale counts/densities at individual cell areas:
    pols, _ = load_polygons(os.path.join(prms.output_dir, 'counts_bin_1.txt'))
    pols_m = convert_to_EPSG(pols, in_epsg=prms.input_epsg, out_epsg=prms.internal_epsg)
    polareas = np.array([pol.area * (prms.epsg_scaling2km ** 2) for pol in pols_m.geoms])  # in km^2
    area_scaling = 1 / prms.density_scaling_factor
    estim.load_densities(inputfile,
                         scaling_factor=area_scaling,
                         rescale_to_polygons_areas=polareas)
    estim.load_bins(prms.bins_file)
    minmags = estim.bins['mins']
    maxmags = estim.bins['maxs']

    # Load (a, b) parameters values:
    # --> NB: values obtained for counts/densities already scaled at individual cell areas!
    if args.use_bootstrap_results is None:
        grt_prms = np.loadtxt(os.path.join(prms.output_dir, 'ab_values.txt'), delimiter=';')
    else:
        grt_prms = np.loadtxt(os.path.join(prms.output_dir, f'ab_values_{args.use_bootstrap_results[0]}.txt'),
                              delimiter=';')

    if estim.cellinfo is None:
        xy = grt_prms[:, :2]  # Longitudes, Latitudes
    else:
        xy = estim.cellinfo[:, :2]
    if prms.fmd_info_file:
        # Load FMD information (Mmin, Mmax, and optionally bin durations):
        estim.load_fmd_info(filename=prms.fmd_info_file, verbose=prms.is_verbose)
    if estim.bin_durations is None:
        estim.bin_durations = np.tile(estim.bins['durations'].reshape((1, estim.nbins)),
                                            (estim.ncells, 1))

    # Plot FMDs for all pixels:
    indices = [i - 1 for i in args.index]  # Transform 1-based indices into 0-based indices
    for indx in indices:
        if prms.bounds_file is not None:
            cellgeom = Point(xy[indx, :].tolist())
            print(f">> Plot FMD for pixel {indx + 1} ({cellgeom.x}, {cellgeom.y})")
        else:
            cellgeom = multipol.geoms[indx]
            print(f">> Plot FMD for polygonal cell {indx + 1}")

        a = _set_value_if_not_nan(grt_prms[indx, 2])
        b = _set_value_if_not_nan(grt_prms[indx, 3])
        da = _set_value_if_not_nan(grt_prms[indx, 4])
        db = _set_value_if_not_nan(grt_prms[indx, 5])
        rho = _set_value_if_not_nan(grt_prms[indx, 6])
        mc = _set_value_if_not_nan(grt_prms[indx, 7])

        if args.no_model:
            # Set a and b values to None to deactivate FMD model plot:
            a = b = None
        else:
            if (a is None) or np.isinf(a):
                a = None
                b = None
                da = None
                print(f'Warning: Invalid a-value detected. Do not plot FMD model.')

        if prms.fmd_info_file:
            mmax = _set_value_if_not_nan(estim.cellinfo[indx, 3], additional_condition=lambda x: x > 0)
        else:
            mmax = None

        # Reconstruct annual seismicity rates for the current cell (normalized to scaling area):
        cell_intensities = np.reshape(estim.densities[indx, 2:], (estim.nbins, ))
        cell_durations = np.reshape(estim.bin_durations[indx, :], (estim.nbins,))

        # Remove unused bins (with durations <= 0):
        ib = (estim.bin_durations[indx, :] > 0)
        if np.all(cell_intensities == 0):
            print(f'Error: Event counts equal to 0 in every magnitude bin. Skipping.')
            continue

        fmd_histogram(minmags[ib],
                      maxmags[ib],
                      cell_intensities[ib] / cell_durations[ib],
                      envelope,
                      cellgeom,
                      a=a,
                      b=b,
                      da=da,
                      modelmc=mc,
                      show_inset=True,
                      mmax=mmax,
                      savefig=os.path.join(args.output_directory, f"fmd_pixel_{indx + 1}.png"),
                      showgrid=args.showgrid)