create_data.py 7.41 KB
Newer Older
felikskiszkurno's avatar
felikskiszkurno committed
1
2
3
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
4
Created on 15.01.2021
felikskiszkurno's avatar
felikskiszkurno committed
5

6
@author: Feliks Kiszkurno
felikskiszkurno's avatar
felikskiszkurno committed
7
8
"""

9
10
11
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
felikskiszkurno's avatar
felikskiszkurno committed
12

13
14
15
16
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.ert as ert

17
18
import slopestabilitytools

19

20
def create_data(test_name, test_config, max_depth):
21
22
23
24
    from main import settings

    world_boundary_v = [-9 * max_depth, 0]  # [NW edge] relatively to the middle
    world_boundary_h = [9 * max_depth, -9 * max_depth]  # [SE edge]
felikskiszkurno's avatar
felikskiszkurno committed
25
26
    # world_boundary_v = [-500, 0]  # [right, left border] relatively to the middle
    # world_boundary_h = [500, -100]  # [top, bottom border]
felikskiszkurno's avatar
felikskiszkurno committed
27
28
29
30
31

    test_results = {}

    # INPUT MODEL - SUBSURFACE START #
    world = mt.createWorld(start=world_boundary_v, end=world_boundary_h,
32
                           layers=test_config['layers_pos'])  # ,
felikskiszkurno's avatar
felikskiszkurno committed
33
34
35
36
37
    # marker=np.linspace(1, tests_horizontal['layer_n']['hor_1'],
    #                  tests_horizontal['layer_n']['hor_1']))

    geometry = world  # +block

felikskiszkurno's avatar
felikskiszkurno committed
38
39
    fig_geometry, ax_geometry = plt.subplots(1)
    pg.show(geometry, ax=ax_geometry)
40
    ax_geometry = slopestabilitytools.set_labels(ax_geometry)
felikskiszkurno's avatar
felikskiszkurno committed
41
    ax_geometry.set_title('1 Geometry of the model')
42
    fig_geometry.tight_layout()
43
44
45
    fig_geometry.savefig('results/figures/png/' + test_name + '_1_geometry.png', bbox_inches="tight")
    #fig_geometry.savefig('results/figures/pdf/' + test_name + '_1_geometry.pdf', bbox_inches="tight")
    #fig_geometry.savefig('results/figures/eps/' + test_name + '_1_geometry.eps', bbox_inches="tight")
felikskiszkurno's avatar
felikskiszkurno committed
46

47
48
49
    measurement_scheme = ert.createERTData(
        elecs=np.linspace(start=-8 * max_depth, stop=8 * max_depth, num=8 * max_depth + 1),
        schemeName='dd')
50

felikskiszkurno's avatar
felikskiszkurno committed
51
52
53
54
55
56
    for electrode in measurement_scheme.sensors():
        geometry.createNode(electrode)
        geometry.createNode(electrode - [0, 0.1])  # What does it do?

    mesh = mt.createMesh(geometry, quality=34)  # , area=2)#

57
    resistivity_map = test_config['rho_values']  # [0]
felikskiszkurno's avatar
felikskiszkurno committed
58

felikskiszkurno's avatar
felikskiszkurno committed
59
60
    fig_model, ax_model = plt.subplots(1)
    pg.show(mesh, data=resistivity_map, label=pg.unit('res'), showMesh=True, ax=ax_model)
61
    ax_model = slopestabilitytools.set_labels(ax_model)
felikskiszkurno's avatar
felikskiszkurno committed
62
    ax_model.set_title('2 Mesh and resistivity distribution')
63
64
    fig_model.tight_layout()
    fig_model.savefig('results/figures/png/' + test_name + '_2_meshdist.png', bbox_inches="tight")
65
66
    #fig_model.savefig('results/figures/pdf/' + test_name + '_2_meshdist.pdf', bbox_inches="tight")
    #fig_model.savefig('results/figures/eps/' + test_name + '_2_meshdist.eps', bbox_inches="tight")
felikskiszkurno's avatar
felikskiszkurno committed
67

felikskiszkurno's avatar
felikskiszkurno committed
68
69
70
71
72
73
74
75
76
77
78
79
80
    input_model = pg.solver.parseMapToCellArray(resistivity_map, mesh)  # rename to input_mesh

    # INPUT MODEL - SUBSURFACE MODEL END ###

    # SIMULATE ERT MEASUREMENT - START ###
    data = ert.simulate(mesh, scheme=measurement_scheme, res=resistivity_map, noiseLevel=1, noiseAbs=1e-6, seed=1337)
    data.remove(data['rhoa'] < 0)
    # SIMULATE ERT MEASUREMENT - END ###

    ert_manager = ert.ERTManager(sr=False, useBert=True, verbose=True, debug=False)

    # RUN INVERSION #
    k0 = pg.physics.ert.createGeometricFactors(data)
81
    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.25, paraMaxCellSize=2, paraDepth=2 * max_depth,
82
83
                                        quality=34, zPower=0.4)

84
    result_full = ert_manager.inv.model
85
86
87
88
89
    # result_array = pg.utils.interperc(result_full, 5)
    # result_lim = result_full.array()
    # result_lim[np.where(result_array > max(resistivity_map[1]))] = float("NaN")
    # result_lim[np.where(result_array < min(resistivity_map[1]))] = float("NaN") # min(resistivity_map[1])
    # result_array = result_lim
90
    result_array = result_full.array()
91
    result_array_norm = slopestabilitytools.normalize(result_array)
92

felikskiszkurno's avatar
felikskiszkurno committed
93
    fig_result, ax_result = plt.subplots(1)
94
95
    pg.show(ert_manager.paraDomain, result_full, label=pg.unit('res'), showMesh=True, ax=ax_result)
    ax_result = slopestabilitytools.set_labels(ax_result)
felikskiszkurno's avatar
felikskiszkurno committed
96
    ax_result.set_title('3 Result')
97
98
    fig_result.tight_layout()
    fig_result.savefig('results/figures/png/' + test_name + '_3_result.png', bbox_inches="tight")
99
100
    #fig_result.savefig('results/figures/pdf/' + test_name + '_3_result.pdf', bbox_inches="tight")
    #fig_result.savefig('results/figures/eps/' + test_name + '_3_result.eps', bbox_inches="tight")
felikskiszkurno's avatar
felikskiszkurno committed
101

102
103
104
    input_model2 = pg.interpolate(srcMesh=mesh, inVec=input_model, destPos=ert_manager.paraDomain.cellCenters())

    input_model2_array = input_model2.array()
105
    input_model2_array_norm = slopestabilitytools.normalize(input_model2_array)
106

felikskiszkurno's avatar
felikskiszkurno committed
107
108
    fig_input, ax_input = plt.subplots(1)
    pg.show(ert_manager.paraDomain, input_model2, label=pg.unit('res'), showMesh=True, ax=ax_input)
109
    ax_input = slopestabilitytools.set_labels(ax_input)
felikskiszkurno's avatar
felikskiszkurno committed
110
    ax_input.set_title('4 Model on inv mesh')
111
112
    fig_input.tight_layout()
    fig_input.savefig('results/figures/png/' + test_name + '_4_modelinv.png', bbox_inches="tight")
113
114
115
116
    #fig_input.savefig('results/figures/pdf/' + test_name + '_4_modelinv.pdf', bbox_inches="tight")
    #fig_input.savefig('results/figures/eps/' + test_name + '_4_modelinv.eps', bbox_inches="tight")

    #if not settings['norm_class']:
felikskiszkurno's avatar
felikskiszkurno committed
117

118
        # Create classes labels
119
120
121
122
123
124
    classes = []
    resistivity_values = []
    for pair in resistivity_map:
        resistivity_values.append(pair[1])
    # print(resistivity_values)

125
    # TODO: This has to be rewritten for more complicated cases
126
127
128
129
130
131
132
133
    for value in input_model2:
        # print(value)
        res_diff = np.abs(value * np.ones_like(resistivity_values) - resistivity_values)
        # print(res_diff)
        res_index = np.argmin(res_diff)
        # print(res_index)
        classes.append(res_index)

134
135
136
137
    #elif settings['norm_class']:

    classesn = slopestabilitytools.assign_classes(input_model2_array_norm)

138
139
140
141
142
143
    # Create sensitivity values
    jac = ert_manager.fop.jacobian()  #
    # Normalization only for visualization!

    # Coverage = cumulative sensitivity = all measurements
    cov = ert_manager.coverage()
144
    # pg.show(ert_manager.paraDomain, cov, label="Logarithm of cumulative sensitivity")
145

146
147
148
149
150
151
152
    rho_arr = []
    for entry in resistivity_map:
        rho_arr.append(entry[1])
    rho_arr = np.array(rho_arr)
    rho_max = np.max(rho_arr)
    rho_min = np.min(rho_arr)

153
    # TODO: this assumes only two resistivities, extend it to consider more
154
155
156
    result_array[np.where(result_array < rho_min)] = rho_min
    result_array[np.where(result_array > rho_max)] = rho_max

157
158
    result_array_norm = slopestabilitytools.normalize(result_array)

159
160
161
162
    experiment_results = pd.DataFrame(data={'X': ert_manager.paraDomain.cellCenters().array()[:, 0],
                                            'Y': ert_manager.paraDomain.cellCenters().array()[:, 1],
                                            'Z': ert_manager.paraDomain.cellCenters().array()[:, 2],
                                            'INM': input_model2_array,
163
                                            'INMN': input_model2_array_norm,
164
                                            'RES': result_array,
165
                                            'RESN': result_array_norm,
166
                                            'SEN': cov,
167
168
                                            'CLASS': classes,
                                            'CLASSN': classesn})
169
170
171
172
173

    # experiment_results.to_csv('results/results/'+test_name+'.csv')

    # test_results[test_name] = experiment_results

174
    experiment_results.to_csv('results/results/' + test_name + '.csv')
175

176
    return experiment_results, rho_max, rho_min