gocad ascii files?

[andrea-bistacchi]

Hello, do you have tools to read and write gocad ascii files, including triangulated surfaces (.ts), polylines (.pl), points (.vs), etc.?

Thanks very much!

[banesullivan]

Can you share documentation (links) for these file types - it likely wouldn't be too difficult to make readers and writers for these formats in PVGeo

[banesullivan]

Looks like this might be the documentation for those file types: http://paulbourke.net/dataformats/gocad/gocad.pdf

[andrea-bistacchi]

Yes, this is it! The structure is quite simple, not very different from an OBJ or PLY. You have: -> A header defining units, CRS, the name of properties, the geological unit and meaning (e.g. fault vs. horizon), etc.

-> Nodes with X Y Z and possibly other properties defined on a per-node basis.

If you have just X Y Z the node lines start with VRTX, then node index, then X Y Z.

VRTX 1 123.456 678.901 234.567 VRTX 2 123.456 678.901 234.567

If there are other properties, you have a PVTX prefix and the other properties follow X Y Z.

PVRTX 1 123.456 678.901 234.567 23 45 PVRTX 2 123.456 678.901 234.567 23 45

-> Connectivity

For polylines:

SEG 1 2 SEG 2 3

For triangulated surfaces:

TRGL 1 2 3 TRGL 2 3 4

For tetrahedral meshes:

TETRA 1 2 3 4 TETRA 2 3 4 5

And obviously there is no connectivity for point clouds (called pointsets in gocad).

-> END - after the END statement you can have another object in the same file. Sometimes you have objects of one kind only (e.g. triangulated surfaces indicated by .ts extension) and some other times the objects can be mixed in the same file.

This is the very basics. The other info in the files are optional (like constraints, borders, border stones = a node at the end of one border), but you might encounter them when importing from some sources.

I can send you some test file if you like. Or we can try writing the reader/writer by ourselves if you have a suitable starting point, e.g. a reader/writer for a similar format.

Thanks very much!

[akaszynski]

If you can provide us with a prototype of the file reader as a pull request, we can start adding to it. pyvista also uses meshio, and you may consider adding it there as it will be able to reach more users that way.

[andrea-bistacchi]

Hello, I am very happy of the interest I have raised for this kind of file format. I think it is justified since the gocad ASCII format is a widely used exchange format between different geomodelling packages (including all the commercial ones).

So I can provide:

1. Fragments of codes in different languages that are able to read gocad-ASCII.

2. A code I developed some time ago to write gocad-ASCII in Matlab. I can convert this to Python, but maybe it is easier for somebody knowing better than me the inner workings of VTK.

3. A sample file including examples of different objects: pointsets, polylines, triangulated surfaces, tetrahedral meshes, and maybe voxets and S-Grids. The latter is a gocad-specific object, given by an hexahedral mesh that is deformed to be conformal to a given stratigraphy. I do not know if this has some counterpart in VTK.

I guess that this is what we need, correct?

[andrea-bistacchi]

I list here all codes I have found online that already do, at least partly, what we need.

ParaViewGeo - an old and apparently no more actively developed spin-off from Paraview, had readers for gocad-ASCII and other formats:

http://paraviewgeo.objectivity.ca/documentation/paraviewgeo-plugins-documentation-1/readers

Coseis has a gocad reader that is supposed to read voxets and tsurfs (triangulated surfaces):

https://github.com/fengw/coseis/blob/master/cst/gocad.py

VtkGocadplot seems to be able to read tsurfs:

http://sepwww.stanford.edu/public/docs/sep111/bob1/paper_html/

SeisSol seems again focused on tsurfs:

https://github.com/SeisSol/Meshing/tree/master/GocadRelatedScripts

Gocad2Flac3d.py reads tetrahedral meshes:

http://www.geo.tu-freiberg.de/~apelm/pub/gocad2flac3d.py

According to me, the simpler approach seems the one implemented here in Coseis:

https://github.com/fengw/coseis/blob/master/cst/gocad.py

[andrea-bistacchi]

Putting together all this, I think we need to parse the gocad-ASCII files to obtain data from three main sections:

-> header:

• name of each object,
• units, color and other cosmetic properties for visualization (line style etc.),
• Coordinate Reference Frame (actually it should be the same for all objects in one file),
• units of Z (time or depth, and meters or feet, seconds or milliseconds),
• positive Z direction (up or down),
• name and dimension of properties exceeding X Y Z (if present - note that properties with dimension 1 are scalars, dimension 3 are vectors and in theory any dimension is allowed).

Then probably there are other fields I do not remember now. The question is: which fields of this kind do we need to be imported in VTK, or for conversion (e.g. positive Z direction up or down).

-> X Y Z coordinates of nodes and other properties recorded on nodes - these lines are marked by VRTX or PVRTX keywords and I guess that they can easily isolated with something like (as in Coseis cited above):

line = lines[counter].strip()
counter += 1
f = line.split()
(...)
elif f[0] in ('VRTX', 'PVRTX'):
vrtx += [[float(f[2]), float(f[3]), float(f[4])]]

-> topology/meshing = connectivity between nodes to obtain cells - these lines are absent in pointsets and are marked as SEG, TRGL, and TETRA for polylines, tsurfs and tetrahedral meshes. I guess that these can be extracted as above.

Finally, we must also look for the END keyword that ends an object and leads to the next one.

Maybe we can start form here for voxets and tsurfs and then expand it for pointsets, polylines, and tetrahedral meshes? I can try but I am sure you know better than me which are supposed to be the target classes where I have to map the data extracted from gocad-ASCII objects (listed above).

[andrea-bistacchi]

Hello I have a prototype gocad reader/writer. Shall I submit it to be included in PyVista or PVGeo?

I do not think Meshio is a good idea since, if we would like to go for a larger library, we should point directly to VTK (more efficient in C++ and there was a gocad reader already in a paraview spin-off called paraviewgeo).

[akaszynski]

I think adding it to pyvista might be best, or you could consider creating a new python package to support it. I'm always wary of feature bloat, and while including it in pyvista would be great, we'll have to maintain it and test it.

Either way, if it's a relatively light library, maybe including it within pyvista won't be a big deal.

[andrea-bistacchi]

Hi and thanks for your comment. I expect to have one function for reading gocad files and one for writing. The gocad ascii format is the same since probably 20 years, so I do not expect that you will have to do a lot of work to “maintain” these functions. However I also understand your point. Maybe I can create a small package, then you can see what is the best way to deal with it.

[banesullivan]

Since gocad is more geoscience oriented, it might make sense to include in PVGeo (which is just an extension package to PyVista at this point).

I'm leaning towards leaving PyVista's supported file types as is (with meshio and adding support for VTK readers that already exist) and having external libraries like PVGeo handle file IO for odd/domain-specific formats.

As for how this might be implemented in PVGeo, @bistek, ping me on slack and I can help you open a PR with the methods in the right place so that we can set up an "Algorithm". Perhaps go ahead and fork PVGeo and add a new directory called gocad with all of the methods you are writing in a new module, then I will help you set up an "algorithm" class that glues it all together so the routines can be used in ParaView and in the typical PVGeo fashion.

Also, do you have files that we can use for testing?

[andrea-bistacchi]

Quick update on this: as suggested by @akaszynski we are carrying out a small independent project on this format. We will release it as soon as it is tested and robust. Best

[GeoMattB]

Hi, this looks like a really interesting tool. Was there a project released? Happy to help with testing.

[RichardScottOZ]

Came across an .sg file today - Structured Grid presumably? Also ascii by the looks.

[RichardScottOZ]

https://github.com/joferkington/python-geoprobe/blob/master/geoprobe/tsurf.py

[RichardScottOZ]

https://programtalk.com/vs2/python/12012/simpeg/SimPEG/Utils/io_utils.py/

[RichardScottOZ]

Having extracted for example, the trianges - can you use them directly?

PVRTX 34384 -10000 6832801.0625 -2500 -99999 -100 -100 -100 -100 PVRTX 34385 -10000 6828283.9375 -2500 -99999 -100 -100 -100 -100 PVRTX 34386 -8506.4375 6830000 -3500 -99999 -100 -100 -100 -100 PVRTX 34387 -11504.375 6830000 -1500 -99999 -100 -100 -100 -100 PVRTX 34388 -10000 6833004.125 -1500 -99999 -100 -100 -100 -100 PVRTX 34389 -10000 6828463.0625 -1500 -99999 -100 -100 -100 -100 PVRTX 34390 -8322.3125 6830000 -2500 -99999 -100 -100 -100 -100 PVRTX 34391 -11139.375 6830000 -500 -99999 -100 -100 -100 -100 PVRTX 34392 -10000 6832818.1875 -500 -99999 -100 -100 -100 -100 PVRTX 34393 -10000 6828814.0625 -500 -99999 -100 -100 -100 -100 PVRTX 34394 -8262.75 6830000 -1500 -99999 -100 -100 -100 -100 PVRTX 34395 -8510.1875 6830000 -500 -99999 -100 -100 -100 -100 TRGL 29713 29716 29715 TRGL 30000 30003 29938 TRGL 19632 19634 19631 TRGL 183 186 185 TRGL 23776 30228 23781 TRGL 342 340 338 TRGL 29996 30056 29999 TRGL 2996 2999 2998 TRGL 1909 1908 1852 TRGL 31344 31342 31343 TRGL 3117 13710 13724 TRGL 22352 23032 23031 TRGL 8170 8171 7189

[RichardScottOZ]

e.g. (1256, 3) (2394, 3) = vertex, triangles <class 'numpy.ndarray'> <class 'numpy.ndarray'> [ 6.85672095e+04 7.04030571e+06 -5.88702338e+03] (mean of vertices) - from https://programtalk.com/vs2/python/12012/simpeg/SimPEG/Utils/io_utils.py/

[RichardScottOZ]

The vertices being fairly obviously PolyData-able.

[RichardScottOZ]

https://energymining.sa.gov.au/minerals/geoscience/geoscientific_data/3d_geological_models/sagrm

SA_Geophysics_Gravity_Inversion_Constrained_IsoSurface_2-63212.zip

[RichardScottOZ]

something like this, from the Slack, where tri-index = triangles faces = np.empty_like(tri_index, shape=(tri_index.shape[0], 4) faces[:, 0] = 3 faces[:, 1:] = tri_index

[RichardScottOZ]

[RichardScottOZ]

The other thing with the above was the faces arrays start at 1, so presumed subtracting one was the thing to do

[RichardScottOZ]

Looking at this again now, clearly! :)

[RichardScottOZ]

@bistek - what are ATOM lines in TS files?

[andrea-bistacchi]

Hello, they should be nodes shared between different objects, used at links or constraints, I guess. Maybe you can skip them.

[RichardScottOZ]

Thanks! Not quite - see here

https://github.com/pyvista/pyvista-support/issues/403

[RichardScottOZ]

https://github.com/RichardScottOZ/GOCAD_TS-Surface-Reader

[RichardScottOZ]

Note mtpy goes back and forth from their model format into sg ascii grids

[RichardScottOZ]

https://github.com/esys-escript/esys-escript.github.io/blob/master/doc/examples/usersguide/voxet_reader.py