create_data.py 6.67 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):
felikskiszkurno's avatar
felikskiszkurno committed
21
22
23
24
    world_boundary_v = [-20 * max_depth, 0]  # [NW edge] relatively to the middle
    world_boundary_h = [20 * max_depth, -5 * max_depth]  # [SE edge]
    # 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
25
26
27
28
29

    test_results = {}

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

    geometry = world  # +block

felikskiszkurno's avatar
felikskiszkurno committed
36
37
    fig_geometry, ax_geometry = plt.subplots(1)
    pg.show(geometry, ax=ax_geometry)
38
    ax_geometry = slopestabilitytools.set_labels(ax_geometry)
felikskiszkurno's avatar
felikskiszkurno committed
39
    ax_geometry.set_title('1 Geometry of the model')
40
41
42
43
    fig_geometry.tight_layout()
    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
44

felikskiszkurno's avatar
felikskiszkurno committed
45
    measurement_scheme = ert.createERTData(elecs=np.linspace(start=-8*max_depth, stop=8*max_depth, num=8*max_depth+1),
46
47
                                           schemeName='dd')

felikskiszkurno's avatar
felikskiszkurno committed
48
49
50
51
52
53
    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)#

54
    resistivity_map = test_config['rho_values']  # [0]
felikskiszkurno's avatar
felikskiszkurno committed
55

felikskiszkurno's avatar
felikskiszkurno committed
56
57
    fig_model, ax_model = plt.subplots(1)
    pg.show(mesh, data=resistivity_map, label=pg.unit('res'), showMesh=True, ax=ax_model)
58
    ax_model = slopestabilitytools.set_labels(ax_model)
felikskiszkurno's avatar
felikskiszkurno committed
59
    ax_model.set_title('2 Mesh and resistivity distribution')
60
61
62
63
    fig_model.tight_layout()
    fig_model.savefig('results/figures/png/' + test_name + '_2_meshdist.png', bbox_inches="tight")
    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
64

felikskiszkurno's avatar
felikskiszkurno committed
65
66
67
68
69
70
71
72
73
74
75
76
77
    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)
felikskiszkurno's avatar
felikskiszkurno committed
78
    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.25, paraMaxCellSize=2, paraDepth=2*max_depth,
79
80
                                        quality=34, zPower=0.4)

81
82
83
84
85
86
87
    result_full = ert_manager.inv.model
    #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
    result_array = result_full.array()
88

felikskiszkurno's avatar
felikskiszkurno committed
89
    fig_result, ax_result = plt.subplots(1)
90
91
    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
92
    ax_result.set_title('3 Result')
93
94
95
96
    fig_result.tight_layout()
    fig_result.savefig('results/figures/png/' + test_name + '_3_result.png', bbox_inches="tight")
    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
97

98
99
100
101
    input_model2 = pg.interpolate(srcMesh=mesh, inVec=input_model, destPos=ert_manager.paraDomain.cellCenters())

    input_model2_array = input_model2.array()

felikskiszkurno's avatar
felikskiszkurno committed
102
103
    fig_input, ax_input = plt.subplots(1)
    pg.show(ert_manager.paraDomain, input_model2, label=pg.unit('res'), showMesh=True, ax=ax_input)
104
    ax_input = slopestabilitytools.set_labels(ax_input)
felikskiszkurno's avatar
felikskiszkurno committed
105
    ax_input.set_title('4 Model on inv mesh')
106
107
108
109
    fig_input.tight_layout()
    fig_input.savefig('results/figures/png/' + test_name + '_4_modelinv.png', bbox_inches="tight")
    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")
felikskiszkurno's avatar
felikskiszkurno committed
110

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
    # Create classes labels
    classes = []
    resistivity_values = []
    for pair in resistivity_map:
        resistivity_values.append(pair[1])
    # print(resistivity_values)

    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)

    # Create sensitivity values
    jac = ert_manager.fop.jacobian()  #
    # Normalization only for visualization!

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

134
135
136
137
138
139
140
141
142
143
    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)

    result_array[np.where(result_array < rho_min)] = rho_min
    result_array[np.where(result_array > rho_max)] = rho_max

144
145
146
147
148
149
150
    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,
                                            'RES': result_array,
                                            'SEN': cov,
                                            'CLASS': classes})
151
152
153
154
155

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

    # test_results[test_name] = experiment_results

156
    experiment_results.to_csv('results/results/' + test_name + '.csv')
157

158
159
160


    return experiment_results, rho_max, rho_min