create_data.py 4.33 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
19
20
def create_data(test_name, test_config, max_depth):
    world_boundary_v = [-10 * max_depth, 0]  # [right, left border] relatively to the middle
    world_boundary_h = [10 * 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
32
33
34
35
36
37
38
    # marker=np.linspace(1, tests_horizontal['layer_n']['hor_1'],
    #                  tests_horizontal['layer_n']['hor_1']))

    geometry = world  # +block

    measurement_scheme = ert.createERTData(elecs=np.linspace(start=-45, stop=45, num=91), schemeName='dd')
    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)#

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

    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)
54
55
    model_inverted = ert_manager.invert(data=data, lam=20, paraDX=0.25, paraMaxCellSize=5, paraDepth=max_depth,
                                        quality=34,
felikskiszkurno's avatar
felikskiszkurno committed
56
                                        zPower=0.4)
felikskiszkurno's avatar
felikskiszkurno committed
57
58
    #result = ert_manager.inv.model
    #result_array = result.array()
felikskiszkurno's avatar
felikskiszkurno committed
59

felikskiszkurno's avatar
felikskiszkurno committed
60
    #input_model2 = pg.interpolate(srcMesh=mesh, inVec=input_model, destPos=ert_manager.paraDomain.cellCenters())
felikskiszkurno's avatar
felikskiszkurno committed
61

felikskiszkurno's avatar
felikskiszkurno committed
62
    ##input_model2_array = input_model2.array()
felikskiszkurno's avatar
felikskiszkurno committed
63

felikskiszkurno's avatar
felikskiszkurno committed
64
65
66
67
68
69
70
    # 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,
    #                                         'INPUT_MODEL': input_model2,
    #                                         'RESULT': result})
71
72
73
74
75
76
77
78
79
80
81
82

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

    # test_results[test_name] = experiment_results

    # Interpolate to grid
    grid = pg.createGrid(x=np.linspace(start=-50, stop=50, num=101),
                         y=-pg.cat([0], pg.utils.grange(0.5, max_depth, n=2 * max_depth + 1)),
                         marker=2)
    input_model3 = pg.interpolate(srcMesh=mesh, inVec=input_model, destPos=grid.cellCenters())
    result_grid = ert_manager.invert(data=data, mesh=grid, lam=20, paraDX=0.25, paraMaxCellSize=5, paraDepth=max_depth, quality=34,
                                    zPower=0.4)
felikskiszkurno's avatar
felikskiszkurno committed
83
84
85
86
87
88
89
90
91
92
93

    class_array = np.ones_like(input_model3) * resistivity_map[-1][0]
    layer_id = 1
    layer_depth_previous = 0

    for depth in test_config['layers_pos']:

        class_array[np.where((grid.cellCenters().array()[:, 1] >= depth) & (
                    grid.cellCenters().array()[:, 1] < layer_depth_previous))] = layer_id
        layer_depth_previous = depth
        layer_id += 1
94
95
96
97
98
99

    experiment_results_grid = pd.DataFrame(data={'X': grid.cellCenters().array()[:, 0],
                                                 'Y': grid.cellCenters().array()[:, 1],
                                                 'Z': grid.cellCenters().array()[:, 2],
                                                 'INM': input_model3.array(),
                                                 'RES': result_grid.array(),
felikskiszkurno's avatar
felikskiszkurno committed
100
                                                 'CLASS': class_array})
101
102
103
104

    experiment_results_grid.to_csv('results/results/' + test_name + '.csv')

    return experiment_results_grid