Skip to content

Commit

Permalink
Merge pull request #679 from DHI/slim-dfsu
Browse files Browse the repository at this point in the history
Dfsu 2.0
  • Loading branch information
ecomodeller authored Mar 25, 2024
2 parents 6cef6af + eb08fb9 commit d905579
Show file tree
Hide file tree
Showing 14 changed files with 537 additions and 443 deletions.
387 changes: 171 additions & 216 deletions mikeio/dfsu/_dfsu.py

Large diffs are not rendered by default.

105 changes: 96 additions & 9 deletions mikeio/dfsu/_layered.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

import numpy as np
from mikecore.DfsuFile import DfsuFile, DfsuFileType
import pandas as pd
from scipy.spatial import cKDTree
from tqdm import trange

Expand All @@ -18,18 +19,102 @@
from .._interpolation import get_idw_interpolant, interp2d
from ..spatial import GeometryFM3D, GeometryFMVerticalProfile, GeometryPoint3D
from ..spatial._FM_utils import _plot_vertical_profile
from ._dfsu import _Dfsu, get_nodes_from_source, get_elements_from_source
from ._dfsu import (
_get_dfsu_info,
get_nodes_from_source,
get_elements_from_source,
_validate_elements_and_geometry_sel,
)

if TYPE_CHECKING:
from ..spatial._FM_geometry_layered import Layer


class DfsuLayered(_Dfsu):
class DfsuLayered:
show_progress = False

def __init__(self, filename: str | Path) -> None:
super().__init__(filename)
info = _get_dfsu_info(filename)
self._filename = info.filename
self._type = info.type
self._deletevalue = info.deletevalue
self._time = info.time
self._timestep = info.timestep
self._geometry = self._read_geometry(self._filename)
# 3d files have a zn item
self._items = self._read_items(self._filename)

def __repr__(self):
out = [f"<mikeio.{self.__class__.__name__}>"]

out.append(f"number of elements: {self.geometry.n_elements}")
out.append(f"number of nodes: {self.geometry.n_nodes}")
out.append(f"projection: {self.geometry.projection_string}")
out.append(f"number of sigma layers: {self.geometry.n_sigma_layers}")
if (
self._type == DfsuFileType.DfsuVerticalProfileSigmaZ
or self._type == DfsuFileType.Dfsu3DSigmaZ
):
out.append(
f"max number of z layers: {self.geometry.n_layers - self.geometry.n_sigma_layers}"
)
if self.n_items < 10:
out.append("items:")
for i, item in enumerate(self.items):
out.append(f" {i}: {item}")
else:
out.append(f"number of items: {self.geometry.n_items}")
if self.n_timesteps == 1:
out.append(f"time: time-invariant file (1 step) at {self.time[0]}")
else:
out.append(
f"time: {str(self.time[0])} - {str(self.time[-1])} ({self.n_timesteps} records)"
)
return str.join("\n", out)

@property
def deletevalue(self) -> float:
"""File delete value"""
return self._deletevalue

@property
def n_items(self) -> int:
"""Number of items"""
return len(self.items)

@property
def items(self) -> list[ItemInfo]:
"""List of items"""
return self._items

@property
def start_time(self) -> pd.Timestamp:
"""File start time"""
return self._time[0]

@property
def n_timesteps(self) -> int:
"""Number of time steps"""
return len(self._time)

@property
def timestep(self) -> float:
"""Time step size in seconds"""
return self._timestep

@property
def end_time(self) -> pd.Timestamp:
"""File end time"""
return self._time[-1]

@property
def time(self) -> pd.DatetimeIndex:
return self._time

@property
def geometry(self) -> GeometryFM3D | GeometryFMVerticalProfile:
return self._geometry

@staticmethod
def _read_items(filename: str) -> list[ItemInfo]:
dfs = DfsuFile.Open(filename)
Expand Down Expand Up @@ -148,7 +233,7 @@ def read(

single_time_selected, time_steps = _valid_timesteps(dfs, time)

self._validate_elements_and_geometry_sel(
_validate_elements_and_geometry_sel(
elements, area=area, layers=layers, x=x, y=y, z=z
)
if elements is None:
Expand All @@ -158,7 +243,7 @@ def read(
or (area is not None)
or (layers is not None)
):
elements = self.geometry.find_index(
elements = self.geometry.find_index( # type: ignore
x=x, y=y, z=z, area=area, layers=layers
)
if len(elements) == 0:
Expand All @@ -179,9 +264,11 @@ def read(
node_ids = geometry.node_ids

item_numbers = _valid_item_numbers(
dfs.ItemInfo, items, ignore_first=self.is_layered
dfs.ItemInfo, items, ignore_first=self.geometry.is_layered
)
items = _get_item_info(
dfs.ItemInfo, item_numbers, ignore_first=self.geometry.is_layered
)
items = _get_item_info(dfs.ItemInfo, item_numbers, ignore_first=self.is_layered)
item_numbers = [it + 1 for it in item_numbers]
if layered_data := hasattr(geometry, "is_layered") and geometry.is_layered:
item_numbers.insert(0, 0)
Expand Down Expand Up @@ -341,8 +428,8 @@ def extract_surface_elevation_from_3d(self, n_nearest: int = 4) -> DataArray:
# make 2d nodes-to-elements interpolator
top_el = self.geometry.top_elements
geom = self.geometry.elements_to_geometry(top_el, node_layers="top")
xye = geom.element_coordinates[:, 0:2]
xyn = geom.node_coordinates[:, 0:2]
xye = geom.element_coordinates[:, 0:2] # type: ignore
xyn = geom.node_coordinates[:, 0:2] # type: ignore
tree2d = cKDTree(xyn)
dist, node_ids = tree2d.query(xye, k=n_nearest)
if n_nearest == 1:
Expand Down
113 changes: 99 additions & 14 deletions mikeio/dfsu/_spectral.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,22 +8,109 @@
from tqdm import trange

from ..dataset import DataArray, Dataset
from ..eum import ItemInfo
from ..dfs._dfs import _get_item_info, _valid_item_numbers, _valid_timesteps
from .._spectral import calc_m0_from_spectrum
from ._dfsu import _Dfsu, get_elements_from_source, get_nodes_from_source
from ._dfsu import (
_get_dfsu_info,
get_elements_from_source,
get_nodes_from_source,
_validate_elements_and_geometry_sel,
)
from ..spatial import (
GeometryFMAreaSpectrum,
GeometryFMLineSpectrum,
GeometryFMPointSpectrum,
)


class DfsuSpectral(_Dfsu):
class DfsuSpectral:
show_progress = False

def __init__(self, filename: str | Path) -> None:
super().__init__(filename)
info = _get_dfsu_info(filename)
self._filename = info.filename
self._type = info.type
self._deletevalue = info.deletevalue
self._time = info.time
self._timestep = info.timestep
self._items = info.items
self._geometry = self._read_geometry(self._filename)

def __repr__(self):
out = [f"<mikeio.{self.__class__.__name__}>"]

if self._type is not DfsuFileType.DfsuSpectral0D:
if self._type is not DfsuFileType.DfsuSpectral1D:
out.append(f"number of elements: {self.geometry.n_elements}")
out.append(f"number of nodes: {self.geometry.n_nodes}")
if self.geometry.is_spectral:
if self.geometry.n_directions > 0:
out.append(f"number of directions: {self.geometry.n_directions}")
if self.geometry.n_frequencies > 0:
out.append(f"number of frequencies: {self.geometry.n_frequencies}")
if self.geometry.projection_string:
out.append(f"projection: {self.geometry.projection_string}")
if self.n_items < 10:
out.append("items:")
for i, item in enumerate(self.items):
out.append(f" {i}: {item}")
else:
out.append(f"number of items: {self.geometry.n_items}")
if self.n_timesteps == 1:
out.append(f"time: time-invariant file (1 step) at {self.time[0]}")
else:
out.append(
f"time: {str(self.time[0])} - {str(self.time[-1])} ({self.n_timesteps} records)"
)
return str.join("\n", out)

@property
def geometry(
self,
) -> GeometryFMPointSpectrum | GeometryFMLineSpectrum | GeometryFMAreaSpectrum:
"""Geometry"""
return self._geometry

@property
def deletevalue(self) -> float:
"""File delete value"""
return self._deletevalue

@property
def n_items(self) -> int:
"""Number of items"""
return len(self.items)

@property
def items(self) -> list[ItemInfo]:
"""List of items"""
return self._items

@property
def start_time(self) -> pd.Timestamp:
"""File start time"""
return self._time[0]

@property
def n_timesteps(self) -> int:
"""Number of time steps"""
return len(self._time)

@property
def timestep(self) -> float:
"""Time step size in seconds"""
return self._timestep

@property
def end_time(self) -> pd.Timestamp:
"""File end time"""
return self._time[-1]

@property
def time(self) -> pd.DatetimeIndex:
return self._time

@staticmethod
def _read_geometry(
filename: str,
Expand Down Expand Up @@ -98,8 +185,8 @@ def _get_spectral_data_shape(
self, n_steps: int, elements: Sized | None, dfsu_type: DfsuFileType
) -> Tuple[Tuple[int, ...], Tuple[int, ...], Tuple[str, ...]]:
dims = [] if n_steps == 1 else ["time"]
n_freq = self.n_frequencies
n_dir = self.n_directions
n_freq = self.geometry.n_frequencies
n_dir = self.geometry.n_directions
shape: Tuple[int, ...] = (n_dir, n_freq)
if n_dir == 0:
shape = (n_freq,)
Expand All @@ -109,21 +196,21 @@ def _get_spectral_data_shape(
read_shape = (n_steps, *shape)
elif dfsu_type == DfsuFileType.DfsuSpectral1D:
# node-based, FE-style
n_nodes = self.n_nodes if elements is None else len(elements)
n_nodes = self.geometry.n_nodes if elements is None else len(elements)
if n_nodes == 1:
read_shape = (n_steps, *shape)
else:
dims.append("node")
read_shape = (n_steps, n_nodes, *shape)
shape = (*shape, self.n_nodes)
shape = (*shape, self.geometry.n_nodes)
else:
n_elems = self.n_elements if elements is None else len(elements)
n_elems = self.geometry.n_elements if elements is None else len(elements)
if n_elems == 1:
read_shape = (n_steps, *shape)
else:
dims.append("element")
read_shape = (n_steps, n_elems, *shape)
shape = (*shape, self.n_elements)
shape = (*shape, self.geometry.n_elements)

if n_dir > 1:
dims.append("direction")
Expand Down Expand Up @@ -202,7 +289,7 @@ def read(
single_time_selected, time_steps = _valid_timesteps(dfs, time)

if self._type == DfsuFileType.DfsuSpectral2D:
self._validate_elements_and_geometry_sel(elements, area=area, x=x, y=y)
_validate_elements_and_geometry_sel(elements, area=area, x=x, y=y)
if elements is None:
elements = self._parse_geometry_sel(area=area, x=x, y=y)
else:
Expand All @@ -214,10 +301,8 @@ def read(

geometry, pts = self._parse_elements_nodes(elements, nodes)

item_numbers = _valid_item_numbers(
dfs.ItemInfo, items, ignore_first=self.is_layered
)
items = _get_item_info(dfs.ItemInfo, item_numbers, ignore_first=self.is_layered)
item_numbers = _valid_item_numbers(dfs.ItemInfo, items)
items = _get_item_info(dfs.ItemInfo, item_numbers)
n_items = len(item_numbers)

deletevalue = self.deletevalue
Expand Down
6 changes: 6 additions & 0 deletions mikeio/spatial/_FM_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -381,6 +381,12 @@ def codes(self, v: np.ndarray) -> None:
raise ValueError(f"codes must have length of nodes ({self.n_nodes})")
self._codes = np.array(v, dtype=np.int32)

@property
def boundary_codes(self):
"""Unique list of boundary codes"""
valid = list(set(self.codes))
return [code for code in valid if code > 0]


class GeometryFM2D(_GeometryFM):
def __init__(
Expand Down
Loading

0 comments on commit d905579

Please sign in to comment.