ERTManager invert is zeroing the z (terrain) data in the input data container - issue

[NoteboomM]

Problem description

I'm loading field ERT data that has elevation data for the electrode positions, which I can check with print(pg.z(data)), and query min, max etc. This survives all my processing (calculating 'k', removing negative rhoa...), but after using the ERTManager invert, print(pg.z(data)) is full of zeroes.

Your environment

--------------------------------------------------------------------------------
  Date: Wed Mar 27 11:49:12 2024 W. Europe Standard Time

                OS : Windows
            CPU(s) : 20
           Machine : AMD64
      Architecture : 64bit
               RAM : 15.6 GiB
       Environment : IPython

  Python 3.10.0 | packaged by conda-forge | (default, Nov 20 2021, 02:18:13)
  [MSC v.1916 64 bit (AMD64)]

           pygimli : 1.5.0
            pgcore : 1.5.0
             numpy : 1.26.4
        matplotlib : 3.8.3
             scipy : 1.12.0
           IPython : 8.18.1
           pyvista : 0.43.4
--------------------------------------------------------------------------------

Steps to reproduce

There's more to the code, but the z data survives. The following sequence reproduces the zeroes. input file: WS_1-6.zip

import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
import os, sys
from pygimli.physics import ert

data = ert.load("WS_1-6.ohm") # Load data to container from BERT format

data['k'] = ert.createGeometricFactors(data, numerical=True)

data['rhoa'] = data['r'] * data['k']

data.remove(data["rhoa"] <= 0) # still some negative values, filtered here

mgr = ert.ERTManager(data) # call the ERTManager method
mgr.invert(verbose=True,
           #paraDX=0.3, paraMaxCellSize=10, paraDepth=20, 
           quality=34,
           lam=200, 
           zWeight=1,
           robust=True,
          )

Expected behavior

It would be useful if this terrain data is preserved, so I can use max and/or min on the terrain information to define y-axis limits for plot(s) of inversion results (the automatic limits are generous with DOI, which is good for QC, but includes a lot of whitespace).

Actual behavior

At the moment, as soon as I run the inversion (even if it fails when I don't remove the negative rhoa), the contents of pg.z changes to an array of zeroes. The showResult output has terrain, so the information is being used to define the mesh, I suppose, but somewhere it appears to wipe the data from the input container.

[NoteboomM]

I'll add a note if it helps with prioritising: There's a simple workaround, simply storing the terrain data separately after ert.load:

data = ert.load("WS_1-6.ohm") # Load data to container from BERT format
terrain = pg.z(data)

Then I can do what I want with that separate array...

[halbmy]

When loading the data, there is no difference made whether the problem is 2D or 3D, in your case the elevation is in the z component of the data. However, if you step into 2D inversion, a 2D mesh needs to be created and for that reason the elevation is moved into the y component of the data where it is preserved and where you should be able to access it. Can you please try and confirm?

[NoteboomM]

Ah, that makes sense...effectively rotating the frame of reference. Yes, pg.y now contains the elevations.

And then the z field would be non-zero for a mesh extruded to 3d, for example?

[halbmy]

Exactly.