Compare datasets

[BenoitURRUTY]

Hi all,

I have some pvtu files and I would like to compare them.

Each file is composed of as a mesh and 15 arrays. I would like to compare an array to another one and keeping the shape of the mesh. Like if I am making a map of the difference.

A dataset looks like this :

+RESULTS:

: UnstructuredGrid (0x7f1c02a6ea68) : N Cells: 1106948 : N Points: 557470 : X Bounds: -2.506e+06, 2.743e+06 : Y Bounds: -2.143e+06, 2.240e+06 : Z Bounds: 0.000e+00, 0.000e+00 : N Arrays: 15

They have the same mesh.

Thanks

[banesullivan]

To clarify, you have one mesh and 15 different data arrays on that mesh which you want to "compare"?

What kinds of comparisons do you want to perform? Each array is accessible by its name using the classic get operator: mesh["array name"]. So you could get all the arrays you want to compare and then perform any array comparison metric you can imagine as the order is the same for each array.

For example, you could have a mesh several arrays and compute the standard deviation between all of the arrays at every node/cell:

import pyvista as pv
from pyvista import examples
import numpy as np

# Loads a sample mesh with data
mesh = examples.download_kitchen()

# Grab specific arrays
names = ["u1", "v1", "w1"]
arrays = np.array([mesh[n] for n in names])

# Compute the comparison metric and add it to the mesh's arrays
mesh["metric"] = np.std(arrays, axis=0)

# Contour and plot
mesh.contour(scalars="metric").plot(scalars="metric")

[BenoitURRUTY]

it's data from a glacier modelling and each dataset contains variable as the velocity.

I would like to compare the velocity from one dataset to the velocity of the other one.

[akaszynski]

Can you zip and upload the dataset?

[BenoitURRUTY]

here you can found one of my dataset : https://filesender.renater.fr/?s=download&token=179cac8c-97c5-4490-823c-7687f6d3463a

my code

import matplotlib.pyplot as plt
import numpy as np
import pyvista as pv
import glob

path= glob.glob('/home/urrutyb/Documents/PhD_Tippacs/real_vtu/run1*1.pvtu')
path=(sorted(path))
data=pv.read(path[0])

Thanks

[akaszynski]

Thanks for the data! So I can demo it, which of the following values do you want to compare?

dict_keys(['groundedmask', 'bedrock', 'zs', 'zb', 'alpha', 'mu', 'pdc_melt', 'h', 'h residual', 'dhdt', 'smb', 'pdc_area', 'ssavelocity', 'lonlat'])

[BenoitURRUTY]

I would like to compare for example 'ssavelocity' with another 'ssavelocity' from another dataset.

[akaszynski]

Since there are many datasets, are you going to compare it based on adjacency, mean, or a single data set to be the "master"?

[akaszynski]

I've built a small example for you:

# path= glob.glob('/home/urrutyb/Documents/PhD_Tippacs/real_vtu/run1*1.pvtu')
# path=(sorted(path))
# data=pv.read(path[0])

import glob
import numpy as np
import pyvista as pv

# load data
files = glob.glob("run1_2/*.vtu")
blocks = pv.MultiBlock([pv.read(f) for f in files])
for i, block in enumerate(blocks):
    block["node_value"] = np.full(block.n_points, i)
combined = blocks.combine()

# plot all data blocks
combined.plot(scalars="node_value", cpos="xy")
              # screenshot='/tmp/blocks.png')

# combined.plot(scalars="ssavelocity", cpos="xy", cmap='Blues')
combined.plot(scalars="ssavelocity", rng=[0, 700], cpos="xy", cmap='Blues')
              # screenshot='/tmp/flow.png')

# plot intensity
pl = pv.Plotter()
pl.add_mesh(blocks[1], scalars='ssavelocity', cmap='Blues')
pl.add_mesh(blocks[10], scalars='ssavelocity', cmap='Blues')
pl.show(screenshot='/tmp/two_blocks.png', cpos='xy')

# plot vectors with mesh in background
# pl = pv.Plotter()
# pl.add_mesh(blocks[1], color='w', opacity=0.2, style='wireframe')
# pl.add_arrows(blocks[1].points, blocks[1].point_arrays['ssavelocity'], mag=20, cmap='Blues')
# pl.add_mesh(blocks[10], color='w', opacity=0.2, style='wireframe')
# pl.add_arrows(blocks[10].points, blocks[10].point_arrays['ssavelocity'], mag=20, cmap='Blues')
# pl.show(cpos='xy')

# plot vectors without mesh
pl = pv.Plotter()
pl.add_arrows(blocks[1].points, blocks[1].point_arrays['ssavelocity'], mag=20, cmap='Blues')
pl.add_arrows(blocks[10].points, blocks[10].point_arrays['ssavelocity'], mag=20, cmap='Blues')
pl.show(cpos='xy')

### compare directions
# these appear to be flow vectors of some sort.  Normalize them so we
# can get a reasonable direction comparison
flow_a = blocks[1].point_arrays['ssavelocity'].copy()
flow_a /= np.linalg.norm(flow_a, axis=1).reshape(-1, 1)
flow_b = blocks[10].point_arrays['ssavelocity'].copy()
flow_b /= np.linalg.norm(flow_b, axis=1).reshape(-1, 1)

# plot normalized vectors
pl = pv.Plotter()
pl.add_arrows(blocks[1].points, flow_a, mag=10000, color='b', label='flow_a')
pl.add_arrows(blocks[10].points, flow_b, mag=10000, color='r', label='flow_b')
pl.add_legend()
pl.camera_position = [(-1044239.3240694795, 354805.0268606294, 484178.24825854995),
                      (-1044239.3240694795, 354805.0268606294, 0.0),
                      (0.0, 1.0, 0.0)]
pl.show()

# flow_a that agrees with the mean flow path of flow_b
agree = flow_a.dot(flow_b.mean(0))
pl = pv.Plotter()
pl.add_mesh(blocks[1], scalars=agree, cmap='bwr', stitle='Flow agreement with block 10')
pl.add_mesh(blocks[10], color='w')
pl.show(cpos='xy')

# compare opposite direction
agree = flow_b.dot(flow_a.mean(0))
pl = pv.Plotter()
pl.add_mesh(blocks[1], color='w')
pl.add_mesh(blocks[10], scalars=agree, cmap='bwr', stitle='Flow agreement with block 10')
pl.show(cpos='xy')

[akaszynski]

Also, since this dataset is cool (dad pun intended), mind if we use it on the pyvista example gallery?

[banesullivan]

Dang, awesome work, @akaszynski! FYI I just put a PR in for enabling log color mapping, which might help with this dataset

[BenoitURRUTY]

No problem to be used in the pyvista example gallery. This modeling is just training for me before to perform more simulations. Maybe just says the data come from an output of http://elmerice.elmerfem.org/

This code looks great will maybe useful.

I think I had some difficulty to explain myself yesterday night. I have a dataset as I sent to you for each year. And I would like to compare one year to the next one.

Here, you can find the next dataset https://filesender.renater.fr/?s=download&token=7f21f220-9f04-d40c-8be2-353387d21d39

My goal will be to compare the run1_2 to the run1_3 for the whole Antarctica but also for each block.

Thanks

[BenoitURRUTY]

I think I found my answer in your code. I just didn't know how to use the "point_arrays" type. pyvista is really great. I was used to paraview but I think pyvista will be more useful. Thanks

# load data
files = glob.glob('/home/urrutyb/Documents/PhD_Tippacs/real_vtu/run1*.vtu')
files=sorted(files)
blocks = pv.MultiBlock([pv.read(f) for f in files])
for i, block in enumerate(blocks):
    block["node_value"] = np.full(block.n_points, i)
combined = blocks.combine()

files1 = glob.glob('/home/urrutyb/Documents/PhD_Tippacs/real_vtu/run2*.vtu')
files1=sorted(files1)
blocks1 = pv.MultiBlock([pv.read(f) for f in files1])
for i, block in enumerate(blocks1):
    block["node_value"] = np.full(block.n_points, i)
combined1 = blocks1.combine()

#make the difference
flow_a = blocks[1].point_arrays['ssavelocity'].copy()
flow_a /= np.linalg.norm(flow_a, axis=1).reshape(-1, 1)
flow_a1 = blocks1[1].point_arrays['ssavelocity'].copy()
flow_a1 /= np.linalg.norm(flow_a1, axis=1).reshape(-1, 1)

flow_diff=flow_a1-flow_a

#plot
pl = pv.Plotter() 
pl.add_mesh(blocks[1], color='w', opacity=0.2, style='wireframe')
pl.add_arrows(blocks[1].points, flow_diff, mag=100, label='flow_diff')
pl.add_legend()
pl.show()

[akaszynski]

Leaving open until we add this to the gallery.