Add PyVista mesh support to Field

[banesullivan]

Resolve #32 (see #32 (comment)) by giving the Field class a way to handle PyVista meshes. After this, I see no advantage for implementing a custom PyVista plotting routine but to just let users do PyVista plotting on their own.

Originally, this function would do everything implicitly without importing the mesh packages... IMO, this isn't a very stable/clean way to implement that as it was checking for attributes on the objects which could be shared across types or the name could be present on an unsupported type that isn't what is expected. I implemented PyVista by bare try/except import and then checking if the passed mesh is a PyVista dataset.

I would HIGHLY recommend implementing type checking for the other meshes and add an alert the user if they are passing an unsupported object.

---

Example usage with PyVista from the code in #36 (comment)

...
grid = pv.UniformGrid()
grid.origin = (329700, 4252600, -2700)
grid.spacing = (250, 250, 50)
grid.dimensions = (60, 75, 100)

# Create the kriging model
krig = krige.Ordinary(fit_model,
                      cond_pos=project["temperature"].points.T, 
                      cond_val=project["temperature"]["temperature (C)"])

krig.mesh(grid, name="temperature (C)", )

...

Another concern. This method does not leverage structured grids... does having a structured grid even help? Because if so, there are a number of PyVista mesh types that can be recast to structured grids which would improve performance

[banesullivan]

This is really messy... we should leverage #60 to handle adding all returned arrays to the PyVista and meshio meshes

[MuellerSeb]

This can now be discussed in #65

[banesullivan]

This is probably the best as is and we can use #65 to change this

[MuellerSeb]

@banesullivan I added a workaround for this problem.

[MuellerSeb]

Please make this PR in the develop branch.

[banesullivan]

the base has been changed to develop

[banesullivan]

If you "squash and merge" on GitHub, this should be good

[banesullivan]

I should probably add some tests to make sure all PyVista meshes work

[LSchueler]

Yeah, I think if you add some tests, we can merge this as is and then, when I finally get along with PR #65, I will automatically be reminded to change this.

[banesullivan]

@LSchueler, will do! And I'll get looking at #65 ASAP too

[banesullivan]

Apologies for doing a rebase and force push... this got outdated quickly.

I'm trying to add tests but having trouble running the tests locally...

[banesullivan]

This same error is happening for like all of the tests:

____________________________________________ TestCondition.test_ordinary ____________________________________________

self = <test_condition.TestCondition testMethod=test_ordinary>

    def test_ordinary(self):
        for Model in self.cov_models:
            model = Model(
>               dim=1, var=0.5, len_scale=2, anis=[0.1, 1], angles=[0.5, 0, 0]
            )

tests/test_condition.py:74:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
gstools/covmodel/base.py:175: in __init__
    self.dim = dim
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <[AttributeError("'Gaussian' object has no attribute '_var'") raised in repr()] Gaussian object at 0x7f9820f3be10>
dim = 1

    @dim.setter
    def dim(self, dim):
        # check if a fixed dimension should be used
        if self.fix_dim() is not None:
            print(self.name + ": using fixed dimension " + str(self.fix_dim()))
            dim = self.fix_dim()
        # set the dimension
        if dim < 1 or dim > 3:
            raise ValueError("Only dimensions of 1 <= d <= 3 are supported.")
        self._dim = int(dim)
        # create fourier transform just once (recreate for dim change)
>       self._sft = SFT(ndim=self.dim, **self.hankel_kw)
E       TypeError: __init__() got an unexpected keyword argument 'alt'

gstools/covmodel/base.py:988: TypeError

any suggestions?

[banesullivan]

Well, I took alt out of the HANKEL_DEFAULT locally in covmodel/base.py and now the tests are running

[banesullivan]

I added a basic test in test_srf.py that uses PyVista PolyData. It is important to note that this method will work for all PyVista mesh types as PyVista's general API is used which is shared across mesh types.

I also noticed the the .mesh() method isn't called anywhere else in the testing suite. This might need to be fully implemented to also test how meshio and ogs5py meshes are handled

[MuellerSeb]

@banesullivan welcome back!! have a look at the new mesh class from #65 .. we plan to include the meshIO stuff there and you are warmly welcome to provide input and wishes there (since we mostly had the vtk format in mind creating this)!

the hankel problem comes in your case from an outdated hankel package (update to 1.0 and it should work!)

maybe we could close this PR in favor of #65 and start over with the new mesh class.

What do you think?

[LSchueler]

@banesullivan Are you fine with closing this PR?

[banesullivan]

All good by me to close this PR! The code will remain here for reference if needed and #59 should take precedence

[MuellerSeb]

@banesullivan I reopened this PR, since we decided, that the mesh class will be a separate class in GSTools, that can be passed to the Field.mesh class like this PR does with pyvista meshes.

I think, the way you proposed here is the best one for now. Sorry for this!

[MuellerSeb]

To solve the problem with multiple returns: We could think about allowing a list of names for the name parameter. We could also care about wrong numbers of names (to also use the standard value of a sole name "field"):

def _names(name, cnt):
    name = [name] if isinstance(name, str) else list(name)[:cnt]
    if len(name) < cnt:
        name += [name[-1] + str(i + 1) for i in range(cnt - len(name))]
    return name

Then we could do something like the following after the field was generated:

fields = [out] if isinstance(out, np.ndarray) else out
for f_name, field in zip(_names(name, len(fields)), fields):
    mesh[f_name] = field

This would result in the names "field" and "field1" when called from a Krige subclass, storing the field and the kriging variance.

[banesullivan]

@MuellerSeb, looks good to me! I like your fix for handling the field names

[MuellerSeb]

@banesullivan thanks for supporting this! 👍