create_data.py 4.81 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
17
import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.ert as ert


18
def create_data(test_name, test_config, max_depth):
felikskiszkurno's avatar
felikskiszkurno committed
19
20
21
22
    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
23
24
25
26
27

    test_results = {}

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

    geometry = world  # +block

felikskiszkurno's avatar
felikskiszkurno committed
34
35
36
37
38
    fig_geometry, ax_geometry = plt.subplots(1)
    pg.show(geometry, ax=ax_geometry)
    ax_geometry.set_title('1 Geometry of the model')
    fig_geometry.savefig('results/figures/png/'+test_name+'_1_geometry.png')

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

felikskiszkurno's avatar
felikskiszkurno committed
42
43
44
45
46
47
    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)#

48
    resistivity_map = test_config['rho_values']  # [0]
felikskiszkurno's avatar
felikskiszkurno committed
49

felikskiszkurno's avatar
felikskiszkurno committed
50
51
52
53
54
    fig_model, ax_model = plt.subplots(1)
    pg.show(mesh, data=resistivity_map, label=pg.unit('res'), showMesh=True, ax=ax_model)
    ax_model.set_title('2 Mesh and resistivity distribution')
    fig_model.savefig('results/figures/png/' + test_name + '_2_meshdist.png')

felikskiszkurno's avatar
felikskiszkurno committed
55
56
57
58
59
60
61
62
63
64
65
66
67
    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
68
    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.25, paraMaxCellSize=2, paraDepth=2*max_depth,
69
70
71
72
73
                                        quality=34, zPower=0.4)

    result = ert_manager.inv.model
    result_array = result.array()

felikskiszkurno's avatar
felikskiszkurno committed
74
75
76
77
78
    fig_result, ax_result = plt.subplots(1)
    pg.show(ert_manager.paraDomain, result, label=pg.unit('res'), showMesh=True, ax=ax_result)
    ax_result.set_title('3 Result')
    fig_result.savefig('results/figures/png/' + test_name + '_3_result.png')

79
80
81
82
    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
83
84
85
86
87
    fig_input, ax_input = plt.subplots(1)
    pg.show(ert_manager.paraDomain, input_model2, label=pg.unit('res'), showMesh=True, ax=ax_input)
    ax_input.set_title('4 Model on inv mesh')
    fig_input.savefig('results/figures/png/' + test_name + '_4_modelinv.png')

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
    # 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")

    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})
118
119
120
121
122

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

    # test_results[test_name] = experiment_results

123
    experiment_results.to_csv('results/results/' + test_name + '.csv')
124

125
    return experiment_results