create_data.py 3.7 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):
19
20
    world_boundary_v = [-20 * max_depth, 0]  # [right, left border] relatively to the middle
    world_boundary_h = [20 * max_depth, -4 * max_depth]  # [top, bottom border]
felikskiszkurno's avatar
felikskiszkurno committed
21
22
23
24
25

    test_results = {}

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

    geometry = world  # +block

32
33
34
    measurement_scheme = ert.createERTData(elecs=np.linspace(start=-5*max_depth, stop=5*max_depth, num=6*max_depth+1),
                                           schemeName='dd')

felikskiszkurno's avatar
felikskiszkurno committed
35
36
37
38
39
40
    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)#

41
    resistivity_map = test_config['rho_values']  # [0]
felikskiszkurno's avatar
felikskiszkurno committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55

    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)
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
    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.1, paraMaxCellSize=2, paraDepth=max_depth,
                                        quality=34, zPower=0.4)

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

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

    input_model2_array = input_model2.array()

    # 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})
96
97
98
99
100

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

    # test_results[test_name] = experiment_results

101
    experiment_results.to_csv('results/results/' + test_name + '.csv')
102

103
    return experiment_results