gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
344 stars 131 forks source link

Problem Export Model ERT Inversion #683

Closed tonapawilliam closed 1 month ago

tonapawilliam commented 1 month ago

Excuse Me Sir, Some time ago, I successfully exported the X Y Res data from the inverting result by using the command

np.savetxt("out1.txt", np.column_stack([pg.x(mgr1.paraDomain.cellCenters()), pg.y(mgr1.paraDomain.cellCenters()), model]))

but recently, I can't use this code anymore with the following error statement

  File ~\miniconda3\envs\pg\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File d:\skripsi\code\perbandingan_qualitymesh hasil inversi.py:87
    np.savetxt("out1.txt", np.column_stack([pg.x(mgr1.paraDomain.cellCenters()), pg.y(mgr1.paraDomain.cellCenters()), model]))

  File ~\miniconda3\envs\pg\lib\site-packages\numpy\lib\shape_base.py:652 in column_stack
    return _nx.concatenate(arrays, 1)

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 7299 and the array at index 2 has size 1

I have tried to update pygimli but there is no change. Is there any way to fix this? Thank You

halbmy commented 1 month ago

The array at index 2, i.e. model, has the length 1. You should be able to find out why. Maybe it is a tuple, or you have a single region? In the latter case, you better use mgr.paraModel(model) which maps back region properties.

tonapawilliam commented 1 month ago

I've tried it but it's weird again, is it possible you can fix my code?

# -*- coding: utf-8 -*-
"""
Created on Mon Mar 25 05:18:47 2024
 Perbedaan ParaMaxCell Size
@author: ASUS
"""
import matplotlib.pyplot as plt
import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
from pygimli.meshtools import polytools as plc
from pygimli.meshtools import quality

"Model 2 Lingkaran"

world = mt.createWorld(start = [0, 0], end = [40, -10], layers = None, marker = 1)
circle1 =  mt.createCircle (pos = [12, -4], radius = 2, marker=2)
circle2 =  mt.createCircle (pos = [28, -4], radius = 2, marker=3)

model =  world + circle1 + circle2

pg.show(model)
#%%
mesh = mt.createMesh(model, quality = 34, smooth =[0,10], area = 0.1, verbose=True)

scheme = ert.createData(elecs = np.linspace(start = 0, stop= 40, num = 41), 
                        schemeName = 'wa' )
rhomap = [[1, 50],
          [2, 10],
          [3, 100]]

data = ert.simulate(mesh = mesh, scheme=scheme, res=rhomap)
data['err'] = ert.estimateError(data)
data.save('data2.dat')

mgr =  ert.ERTManager(data, verbose = True)
inv = mgr.invert(data, lam =10, paraDepth = 10, paraMaxCellSize = 0.1, zWeight = 0.02 )

mgr.showResult ()

#%%
np.savetxt("PRM 01.dat", np.column_stack([pg.x(mgr.paraDomain.cellCenters()), 
                                                      pg.y(mgr.paraDomain.cellCenters()), mgr.paraModel(model)]))
Traceback (most recent call last):

  File ~\miniconda3\envs\pg\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File d:\skripsi\code\untitled2.py:45
    pg.y(mgr.paraDomain.cellCenters()), mgr.paraModel(model)]))

  File ~\miniconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py:726 in paraModel
    return self.fop.paraModel(model)

  File ~\miniconda3\envs\pg\lib\site-packages\pygimli\frameworks\modelling.py:635 in paraModel
    mod = model[self.paraDomain.cellMarkers()]

  File ~\miniconda3\envs\pg\lib\site-packages\pygimli\core\mesh.py:138 in __Mesh_getVal
    if self.haveData(key):

ArgumentError: Python argument types in
    None.haveData(Mesh, IVector)
did not match C++ signature:
    haveData(GIMLI::Mesh {lvalue}, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > name)
tonapawilliam commented 1 month ago

I am sorry Sir, I got the problem, I change the line with this code

out = np.column_stack((mgr.paraDomain.cellCenters(), mgr.model))
np.savetxt("Outt.txt", out)

or I change the "model" to mgr.model just like this np.savetxt("out1.txt", np.column_stack([pg.x(mgr.paraDomain.cellCenters()), pg.y(mgr.paraDomain.cellCenters()), mgr.model]))

halbmy commented 1 month ago

Correct. In the original code, model is (other than in many examples) a PLC and not the inversion result (model=mgr.invert()). Nevertheless, the inversion result is story in mgr.model.