search

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
136 stars 53 forks source link

cannot be more sigma than totlayers #605

Closed otzi5300 closed 7 months ago

otzi5300 commented 9 months ago

File read by mikecore can sometimes contain more sigma than total number of layers.

ecomodeller commented 9 months ago

@otzi5300 could you explain the background behind this PR in a bit more detail, and include the output section from the pfs?

It might be useful information for e.g. @JesperGr, since this seems like a bug in the FM engine.

otzi5300 commented 9 months ago

It seems like MIKE is writing a dfsu file with the wrong numbers of sigma layers when selecting 3d output with single layer selected. DFSU file in question have 3 sigma and 24 z layers and as output is 3D field as volume series with layer range 27-27 (top layer).

df = mikeio.read("../oil2d.dfsu") print(df)

dims: (time:7, element:27597) time: 2023-08-24 00:00:00 - 2023-08-24 06:00:00 (7 records) geometry: Flexible Mesh Geometry: Dfsu3DSigmaZ number of nodes: 34104 number of elements: 27597 number of layers: 1 number of sigma layers: 3 Part of output section in question (in the m3fm file): ``` [OUTPUT_3] Touched = 1 include = 1 title = '2d' file_name = 'oil2d.dfsu' type = 2 format = 3 flood_and_dry = 2 coordinate_type = 'PROJCS["SWEREF99_TM",GEOGCS["GCS_SWEREF99",DATUM["D_SWEREF99",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]' zone = 0 input_file_name = || input_format = 1 interpolation_type = 1 use_end_time = 1 first_time_step = 0 last_time_step = 960 time_step_frequency = 120 first_particle = 1 last_particle = 0 particle_frequency = 1 number_of_points = 1 [POINT_1] name = 'Point 1' x = 616915.2324015191 y = 6599389.59576897 z = -30.1403465270996 layer = 1 EndSect // POINT_1 ... [AREA] number_of_points = 4 [POINT_1] x = 558034.9299406767 y = 6566761.336970681 EndSect // POINT_1 [POINT_2] x = 558034.9299406767 y = 6632017.854567259 EndSect // POINT_2 [POINT_3] x = 675795.5348623614 y = 6632017.854567259 EndSect // POINT_3 [POINT_4] x = 675795.5348623614 y = 6566761.336970681 EndSect // POINT_4 layer_min = 27 layer_max = 27 orientation = 0.0 x_origo = 559189.445675203 x_ds = 6076.398602770108 x_npoints = 20 y_origo = 6567401.10675104 y_ds = 6076.398602770108 y_npoints = 12 z_origo = -60.2806930541992 z_ds = 8.037425740559891 z_npoints = 10 EndSect // AREA ```