oceanmodeling / xarray-selafin

An xarray engine for opening Selafin files (TELEMAC)
https://pypi.org/project/xarray-selafin/
The Unlicense
4 stars 2 forks source link

support for 1D result files #27

Closed tomsail closed 8 months ago

tomsail commented 8 months ago

This concerns the particular case of 1D time series generated with the keywords:

TIME SERIES FILE 1             = 'results_1D.slf'
TIME SERIES COORDINATES FILE 1 = 'station.in'

in T2D. (I haven't tested yet in T3D or TOMAWAC)

I can provide the file generated: results_1D.zip (to unzip)

opening it with PyTelTools Serafin gives me the following traceback:

ds = xr.open_dataset('results_1D.slf', engine = 'selafin')
SERAFIN VALIDATION ERROR: Unknown mesh type
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/tomsail/work/python/pyPoseidon/.venv/lib/python3.10/site-packages/xarray/backends/api.py", line 573, in open_dataset
    backend_ds = backend.open_dataset(
  File "/home/tomsail/work/python/pyPoseidon/.venv/lib/python3.10/site-packages/xarray_selafin/xarray_backend.py", line 268, in open_dataset
    slf = read_serafin(filename_or_obj, lang)
  File "/home/tomsail/work/python/pyPoseidon/.venv/lib/python3.10/site-packages/xarray_selafin/xarray_backend.py", line 22, in read_serafin
    resin.read_header()
  File "/home/tomsail/work/python/pyPoseidon/.venv/lib/python3.10/site-packages/xarray_selafin/Serafin.py", line 1020, in read_header
    self.header.from_file(self.file, self.file_size)
  File "/home/tomsail/work/python/pyPoseidon/.venv/lib/python3.10/site-packages/xarray_selafin/Serafin.py", line 885, in from_file
    self._check_dim()
  File "/home/tomsail/work/python/pyPoseidon/.venv/lib/python3.10/site-packages/xarray_selafin/Serafin.py", line 198, in _check_dim
    raise SerafinValidationError("Unknown mesh type")
xarray_selafin.Serafin.SerafinValidationError: Unknown mesh type

neither BlueKenue can open it.

Although opening it with the Selafin class works without problem :

slf= Selafin('results_1D.slf')
>>> slf.meshx
array([13.55799961, 13.23700047, 12.07441902, ...,  9.87355423,
        5.1170001 ,  5.99513721])
>>> slf.meshy
array([-12.33300018,  -8.7869997 , -15.15768051, ...,  57.59385681,
        61.93299866,  53.44028473])
>>> slf.get_values(0)
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)
>>> slf.ikle3
array([[   0],
       [   1],
       [   2],
       ...,
       [1146],
       [1147],
       [1148]])
>>> slf.iparam
array([   1,    0,    0,    0,    0,    0,    0, 1149,    0,    1])

The problem is that the function

def _check_dim(self):
        # verify data consistence and determine 2D or 3D
        if self.is_2d:
            if self.nb_nodes_per_elem != 3:
                raise SerafinValidationError("Unknown mesh type")
        else:
            if self.nb_nodes_per_elem != 6:
                raise SerafinValidationError(
                    "The number of nodes per element is not equal to 6"
                )
            if self.nb_planes < 2:
                raise SerafinValidationError("The number of planes is less than 2")

does not take this case in consideration. While the header class in Selafin might be more forgiving?

This might also be relevant to add this issue directly in PyTelTools.

I let you decide @lucduron

lucduron commented 8 months ago

The fix for 2D seems simple, I can provide it. Can you provide a file in 3D with values at elements (if it is used)?