CODEWITHSUNDEEP
numpyparaview
user2398300
user2398300

Reputation: 25

Visualize a 3D numpy vector field in Paraview

I wrote some code that calculates electric and magnetic vector fields for variety of charge, current and geometric configurations. The output is saved in csv format using the numpy savetxt method:

    fieldMap = np.stack((X.ravel(), Y.ravel(), Z.ravel(), Bx.ravel(), By.ravel(), Bz.ravel(), Ex.ravel(), Ey.ravel(), Ez.ravel()))   
    heading = f" x y z Bx By Bz Ex Ey Ez \n"
    np.savetxt(path , fieldMap,  header = heading , comments="" )

With the expected output

x y z Bx By Bz Ex Ey Ez 
3.840000000000000357e-02 3.840000000000000357e-02 7.439999999999999392e-02 3.731101704388582415e-04 3.731101704388588378e-04 -1.787606656866691209e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -0.000000000000000000e+00
3.840000000000000357e-02 3.840000000000000357e-02 5.924444444444443986e-02 2.135196586408618616e-04 2.135196586408621869e-04 -4.585864781494742551e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -4.034679354447678890e-01
3.840000000000000357e-02 3.840000000000000357e-02 4.408888888888888580e-02 -1.031604762420866280e-04 -1.031604762420863569e-04 -4.494106267203430705e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -1.180245238916466377e+00
3.840000000000000357e-02 3.840000000000000357e-02 2.893333333333333174e-02 -2.527944858538051917e-04 -2.527944858538058964e-04 -1.402243480967241794e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -1.861414658185234217e+00
3.840000000000000357e-02 3.840000000000000357e-02 1.377777777777777768e-02 -1.654374219591531604e-04 -1.654374219591537025e-04 1.197934198869413612e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -2.366509916060799412e+00

This format is required for input to a Geant4 particle transport model, and it works as expected. It also plots out as expected in matplotlib. But trying to get the data into Paraview for a better visualization has been remarkably difficult. I first tried just importing the csv file directly but get the following error:

ERROR: In ./VTK/Common/DataModel/vtkTable.cxx, line 355
vtkTable (0x6309affe5ef0): Column "" must have 1000 rows, but has 0.

If I continue and try to run a filter, the program just crashes. It seems that some invisible column with zero rows is being being inserted in the data. Here's a screen shot of the csv table in paraview.

enter image description here

After that failed, I tried conversion to a vtk file format using: pyevtk, the vtk.util.numpy_support Module, converting to vtkImagedata, numpy points array to vtkPolydata, and PyVista. In every case I get the following import error in Paraview.

ERROR: In ./VTK/IO/XMLParser/vtkXMLParser.cxx, line 375
vtkXMLDataParser (0x65093f5e52a0): Error parsing XML in stream at line 19, column 36, byte index 1576: not well-formed (invalid token)

ERROR: In ./VTK/IO/XML/vtkXMLReader.cxx, line 521
vtkXMLImageDataReader (0x65094191b3a0): Error parsing input file.  ReadXMLInformation aborting.

the offending line, in this case 19, is always the second line of the binary block, after the tag in the XML. I finally ran one of the examples in the pyevtk source listing and got the same error, so its not specific to my code.

My first question, is there something less hostile than Paraview that I can use? I just need to run a few visualizations. I'm astounded that something so fundamental as viewing a vector field should be this difficult. I've been working this on and off for over a month but now its the last obstacle to finishing my project. I started coding in FORTRAN on punch cards in 1980 and this is the hardest stop I’ve ever had in programming. Any help would be appreciated.

I already listed everything I've tried. As for expectations, I expected that it would be relatively straight forward to move data between two such widely used software platforms. As for versions, I've tried paraview 5.10.0-RC1 and 5.12.0, Python 3.10 and 3.11.4, numpy 1.21.5 and 1.26.4. I didn't include all the code i tried because the results are so consistent. Any attempt at converting to a vtk format frsults in a xml parsing error, and importing as csv crashes the program everytime I try to filter the table data.

Upvotes: 1

Views: 175

Answers (1)

C-3PO
C-3PO

Reputation: 1213

I encountered a similar challenge a few days ago. The best solution I found was to save the vector-data as a tecplot file, and import it to Paraview.

This is the full code, including the data given in the example:

original_data = '''
x y z Bx By Bz Ex Ey Ez 
3.840000000000000357e-02 3.840000000000000357e-02 7.439999999999999392e-02 3.731101704388582415e-04 3.731101704388588378e-04 -1.787606656866691209e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -0.000000000000000000e+00
3.840000000000000357e-02 3.840000000000000357e-02 5.924444444444443986e-02 2.135196586408618616e-04 2.135196586408621869e-04 -4.585864781494742551e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -4.034679354447678890e-01
3.840000000000000357e-02 3.840000000000000357e-02 4.408888888888888580e-02 -1.031604762420866280e-04 -1.031604762420863569e-04 -4.494106267203430705e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -1.180245238916466377e+00
3.840000000000000357e-02 3.840000000000000357e-02 2.893333333333333174e-02 -2.527944858538051917e-04 -2.527944858538058964e-04 -1.402243480967241794e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -1.861414658185234217e+00
3.840000000000000357e-02 3.840000000000000357e-02 1.377777777777777768e-02 -1.654374219591531604e-04 -1.654374219591537025e-04 1.197934198869413612e-04 -0.000000000000000000e+00 -0.000000000000000000e+00 -2.366509916060799412e+00
'''

original_data  =  [x for x in original_data.split('\n') if x]
header         =  original_data[0].split()
rows           =  [row.split() for row in original_data[1:]]

filename       =  'vector_field_electrical.tec'

with open(filename,'w') as f:
    f.write('VARIABLES = "' + '", "'.join(header) + '"' + '\n')
    f.write(f"ZONE I={len(rows)}, DATAPACKING=BLOCK"    + '\n')
    for j in range(len(rows[0])):
        f.write(" ".join([str(row[j]) for row in rows]) + '\n')

The tecplot file is immediately recognized by Paraview, and everything is imported:

enter image description here

PS. Tecplot writer as a function:

# ---------------- tecplot writer ----------------
def write_tec3d_line(filename, data):
    assert [x.lower() for x in list(data.keys())[:3]] == ['x', 'y', 'z']
    num_points  =  len(list(data.values())[0])
    with open(filename,'w') as f:
        f.write('VARIABLES = "' + '", "'.join(list(data.keys())) + '"' + '\n')
        f.write(f"ZONE I={num_points}, DATAPACKING=BLOCK"    + '\n')
        for column in data.values():
            f.write(" ".join([str(val) for val in column]) + '\n')

# ---------------- main runner ----------------
def main_runner():
    # ---------------- data definition ----------------
    import numpy as np
    X,Y,Z,Bx,By,Bz,Ex,Ey,Ez = list(map(np.array,[[0.038400000000000004, 0.038400000000000004, 0.038400000000000004, 0.038400000000000004, 0.038400000000000004],
                                                [0.038400000000000004, 0.038400000000000004, 0.038400000000000004, 0.038400000000000004, 0.038400000000000004],
                                                [0.0744, 0.05924444444444444, 0.044088888888888886, 0.028933333333333332, 0.013777777777777778], 
                                                [0.00037311017043885824, 0.00021351965864086186, -0.00010316047624208663, -0.0002527944858538052, -0.00016543742195915316], 
                                                [0.00037311017043885884, 0.0002135196586408622, -0.00010316047624208636, -0.0002527944858538059, -0.0001654374219591537], 
                                                [-0.00017876066568666912, -0.00045858647814947426, -0.00044941062672034307, -0.00014022434809672418, 0.00011979341988694136], 
                                                [-0.0, -0.0, -0.0, -0.0, -0.0], [-0.0, -0.0, -0.0, -0.0, -0.0], 
                                                [-0.0, -0.4034679354447679, -1.1802452389164664, -1.8614146581852342, -2.3665099160607994]]))

    # ---------------- write output ----------------
    filename  =  'vector_field_electrical.tec'
    data      =  dict(x  = X.ravel()  , # x,y,z must be the first entries of the dictionary
                      y  = Y.ravel()  , 
                      z  = Z.ravel()  , 
                      Bx = Bx.ravel() , 
                      By = By.ravel() , 
                      Bz = Bz.ravel() , 
                      Ex = Ex.ravel() , 
                      Ey = Ey.ravel() , 
                      Ez = Ez.ravel() )
    write_tec3d_line(filename, data)

# ---------------- main call ----------------
if __name__ == '__main__':
    main_runner()

Upvotes: 0

Related Questions