DHI / mikeio

Read, write and manipulate dfs0, dfs1, dfs2, dfs3, dfsu and mesh files.
https://dhi.github.io/mikeio
BSD 3-Clause "New" or "Revised" License
139 stars 57 forks source link

Faulty conversion of direction on read() in _spectral.py #724

Closed buchmann-dhi closed 2 months ago

buchmann-dhi commented 2 months ago

_spectral.py read() converts rad2deg In _spectral.py, directions are automatically (and always) converted from radian to degree with:

directions = None if dir is None else dir * (180 / np.pi)

https://github.com/DHI/mikeio/blob/7a01cc74d92470ba6ecc1d6b86fb94c478ba7051/mikeio/dfsu/_spectral.py#L138

However, sometimes (depending on Mike user settings), output may already be in degrees. In this case, the conversion is unwanted resulting in much too large values for the directions.

It seems necessary to test the unit type before blindly applying rad2deg.

Reproduction:

fil=r'BC_Globalv3SW_GulfBoundE_2024090512F000.dfsu'
dfs = mikeio.open(fil)
print(dfs.geometry.directions)

Expected output:

[    0.       572.9578  1145.9156  1718.8734  2291.8313  2864.789
  3437.7468  4010.7046  4583.6626  5156.62    5729.578   6302.5356
  6875.4937  7448.4517  8021.409   8594.367   9167.325   9740.282
 10313.24   10886.198  11459.156  12032.114  12605.071  13178.03
 13750.987  14323.945  14896.903  15469.86   16042.818  16615.775
 17188.734  17761.691  18334.65   18907.607  19480.564  20053.523 ]

System information:

jsmariegaard commented 2 months ago

Thanks @buchmann-dhi

After discussion with @ecomodeller, I suggest a solution like this:

Just before the above mentioned line...

This could easily be encapsulated in a function _is_directions_in_rad() for better readability.

Notes: int(EUMUnit.degree) = 2401 int(EUMUnit.radian) = 2400

buchmann-dhi commented 2 months ago

Sounds good @jsmariegaard I would maybe suggest a small change to the suggestion: Define a conversion factor to be used in any case. If the input is already degree (2401), then set it to unity. If input is 2400, then use 180/np.pi. This opens for the possibility of other conversions (although rad and deg are probably all that will be needed). Comments?