Update to use recent MIKE IO and ModelSkill

.github/workflows/full_test.yml

@@ -15,8 +15,8 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-latest, windows-latest]
-        python-version: [3.8, 3.9, "3.10"]
+        os: [ubuntu-latest]
+        python-version: ["3.10","3.13"]
       fail-fast: false
 
     steps:

.github/workflows/notebooks_test.yml

@@ -16,7 +16,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v2
       with:
-        python-version: '3.9'    
+        python-version: '3.13'    
     - name: Install fmdap
       run: |
         pip install .[test,notebooks]    

.github/workflows/scheduled_test.yml

@@ -13,8 +13,8 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-latest, windows-latest]
-        python-version: [3.8, 3.9, "3.10"]
+        os: [ubuntu-latest]
+        python-version: ["3.13"]
 
     steps:
     - uses: actions/checkout@v2

.gitignore

@@ -129,4 +129,4 @@ dmypy.json
 .pyre/
 
 /data/
-/tmp/
+.envrc

docs/index.rst

@@ -1,8 +1,3 @@
-.. fmskill documentation master file, created by
-   sphinx-quickstart on Thu Mar 25 08:11:16 2021.
-   You can adapt this file completely to your liking, but it should at least
-   contain the root `toctree` directive.
-
 .. image:: https://raw.githubusercontent.com/DHI/fmdap/main/images/logo/mike-fm-dap-rgb.svg
 
 |

docs/vision.rst

@@ -7,7 +7,7 @@ Vision
 
 Design principles
 -----------------
-Like `MIKE IO <https://github.com/DHI/mikeio>`_ and `FMSkill <https://github.com/DHI/fmskill>`_, FMDAp should be: 
+Like `MIKE IO <https://github.com/DHI/mikeio>`_ and `ModelSkill <https://github.com/DHI/modelskill>`_, FMDAp should be: 
 
 * Easy to use
 * Easy to install

fmdap/AR1.py

@@ -4,7 +4,7 @@
 from statsmodels.tsa.arima_process import ArmaProcess
 
 
-def phi_to_halflife(phi, dt=1):
+def phi_to_halflife(phi, dt=1.0):
     """Convert the AR(1) propagation parameter phi to half-life"""
     rho = dt / (np.log2(1 / phi))
     return rho
@@ -37,10 +37,10 @@ def estimate_AR1_halflife(df):
     phi = res.params[1]
 
     try:
-        dt = df.index.freq.delta.seconds
+        dt = pd.Timedelta(df.index.freq).total_seconds()
         # dt = pd.infer_freq(df.index).delta.seconds
     except:
-        dt = 1
+        dt = 1.0
     return phi_to_halflife(phi, dt=dt)
 
 

fmdap/diagnostic_collection.py

@@ -3,7 +3,7 @@
 import numpy as np
 import warnings
 import pandas as pd
-from fmskill.comparison import ComparerCollection
+import mikeio
 import fmdap.diagnostic_output
 from fmdap import pfs_helper as pfs
 
@@ -29,7 +29,7 @@ def __init__(self, diagnostics=None, names=None, attrs=None):
 
     @classmethod
     def from_pfs(cls, pfs_file, folder=None, types=[1, 2]):
-        df, DA_type = cls._parse_pfs(pfs_file, types)
+        df, DA_type = cls._parse_PfsDocument(pfs_file, types)
         df = cls._check_file_existance(df, folder)
         dc = cls()
         for _, row in df.iterrows():
@@ -41,16 +41,10 @@ def from_pfs(cls, pfs_file, folder=None, types=[1, 2]):
         return dc
 
     @classmethod
-    def _parse_pfs(cls, pfs_file, types=[1, 2]):
-
-        warnings.filterwarnings("ignore", message="Support for PFS files")
+    def _parse_PfsDocument(cls, pfs_file, types=[1, 2]):
         assert os.path.exists(pfs_file)
-        d = pfs.pfs2dict(pfs_file).get("DATA_ASSIMILATION_MODULE")
-        if d is None:
-            raise ValueError(
-                "'DATA_ASSIMILATION_MODULE' section could not be found in pfs file!"
-            )
-        DA_type = d.get("METHOD", {}).get("type", 0)
+        d = mikeio.read_pfs(pfs_file).DA
+        DA_type = d.METHOD.type
 
         dfd = pfs.get_diagnostics_df(d)
 
@@ -187,7 +181,7 @@ def _get_diag_id(self, diag):
                 diag_id = diag
             else:
                 raise IndexError(
-                    f"diagnostic id {diag} is out of range (0, {len(self)-1})"
+                    f"diagnostic id {diag} is out of range (0, {len(self) - 1})"
                 )
         else:
             raise KeyError("must be str or int")
@@ -299,7 +293,7 @@ def _diagnostics_attribute_to_frame(self, attr):
 
     def skill(self, **kwargs):
         s = self.comparer.skill(**kwargs)
-        s.df = self._split_skill_index(s.df)
+        s.data = self._split_skill_index(s.data)
         return s
 
     def _split_skill_index(self, df):
@@ -317,7 +311,7 @@ def _split_skill_index(self, df):
             return df.set_index(["observation", "selection"])
 
     def scatter(self, **kwargs):
-        return self.comparer.scatter(**kwargs)
+        return self.comparer.plot.scatter(**kwargs)
 
     @property
     def comparer(self):
@@ -326,14 +320,18 @@ def comparer(self):
         return self._comparer
 
     def _get_comparer(self):
-        cc = ComparerCollection()
+        cmps = []
         for n in self.names:
             try:
                 diag = self.diagnostics[n]
-                cc.add_comparer(diag.comparer)
+                cmps.append(diag.comparer)
             except AttributeError:
                 pass
                 # warnings.warn(f"Could not add 'comparer' from {n}. No such attribute.")
+        # TODO this can probably be done cleaner
+        cc = cmps[0]
+        for cmp in cmps[1:]:
+            cc = cc + cmp
         return cc
 
     def iplot(self, title=None, **kwargs):  # pragma: no cover
@@ -383,4 +381,3 @@ def _iplot_add_subplot(fig, row, diag):  # pragma: no cover
             fig.add_hline(y=0.0)
 
         fig["layout"][f"yaxis{row}"]["title"] = diag.eumText
-

fmdap/diagnostic_output.py

@@ -649,7 +649,7 @@ def skill(self, **kwargs):
         return self.comparer.skill(**kwargs)
 
     def scatter(self, **kwargs):
-        return self.comparer.scatter(**kwargs)
+        return self.comparer.plot.scatter(**kwargs)
 
     @property
     def comparer(self):
@@ -658,24 +658,20 @@ def comparer(self):
         return self._comparer
 
     def _get_comparer(self):
-        import fmskill
+        import modelskill as ms
 
         if self.is_point:
-            mod = fmskill.ModelResult(self.df[["Mean_State"]])
-            obs = fmskill.PointObservation(self.df[["Measurement"]], name=self.name)
+            mod = ms.model_result(self.df[["Mean_State"]])
+            obs = ms.PointObservation(self.df[["Measurement"]], name=self.name)
         else:
             df = self.df.reset_index(["x", "y"])
-            mod = fmskill.ModelResult(df[["x", "y", "Mean_State"]], type="track")
-            obs = fmskill.TrackObservation(
+            mod = ms.TrackModelResult(df[["x", "y", "Mean_State"]])
+            obs = ms.TrackObservation(
                 df[["x", "y", "Measurement"]], name=self.name
             )
 
-        con = fmskill.Connector(obs, mod)
-        cc = con.extract()
-        if (cc is None) or (len(cc) == 0):
-            return cc  # None
-        else:
-            return cc[0]
+        return ms.match(obs,mod)
+
 
 
 class MeasurementPointDiagnostic(_DiagnosticIndexMixin, MeasurementDiagnostic):

fmdap/diagnostic_output_altimetry.py

@@ -7,18 +7,7 @@
 
 
 class DiagnosticOutputAltimetry:
-    dfd = None  # dataframe of diagnostic data
-    dfo = None  # dataframe of observations
-    dfda = None  # dataframe of DA steps
-    dfqa = None  # dataframe of non-DA steps
-    is_DA = None
-    nc1 = None  # boundary nodes
-    msh = None  # model mesh
-
-    def __init__(self):
-        self._dfs = mikeio.Dfs0()
-
-    def read(self, file_diag, file_obs, obs_col_name="adt_dhi"):
+    def __init__(self, file_diag, file_obs, obs_col_name="adt_dhi"):
         """Read diagnostic output dfs0 and associated observation dfs0 and store as data frames
 
         Arguments:
@@ -67,7 +56,7 @@ def remove_points_outside_mesh(self):
             raise ValueError("mesh has not been provided. Please use read_mesh() first")
 
         xy = np.vstack([self.dfo.iloc[:, 0].values, self.dfo.iloc[:, 1].values]).T
-        inside = self.msh.contains(xy)
+        inside = self.msh.geometry.contains(xy)
         self.dfo = self.dfo[inside]
 
     def process(self, col_no_DA="mike_wl", track_split_time=3):

fmdap/pfs.py

@@ -1,24 +1,20 @@
-# from collections import namedtuple
-import warnings
+from pathlib import Path
 import numpy as np
 import pandas as pd
 import mikeio
 
 
 class Pfs:
-    def __init__(self, pfs_file=None) -> None:
+    def __init__(self, pfs_file: str | Path) -> None:
         self.d = None
         self._sections = None
         self._model_errors = None
         self._measurements = None
         self._diagnostics = None
 
         if pfs_file:
-            with warnings.catch_warnings():
-                warnings.simplefilter("ignore")
-                pfs = mikeio.Pfs(pfs_file)
-            self.data = pfs.data  # PfsSection
-            self.d = pfs.data.to_dict()  # dictionary
+            pfs = mikeio.PfsDocument(pfs_file)
+            self.d = pfs.targets[0].to_dict()
 
     @property
     def dda(self):
@@ -57,7 +53,7 @@ def measurement_positions(self):
     def validate_positions(cls, mesh, df):
         """Determine if positions are inside mesh and find nearest cell centers"""
         # TODO: handle empty positions
-        assert isinstance(mesh, (mikeio.Mesh, mikeio.dfsu._Dfsu))
+        assert isinstance(mesh, (mikeio.Mesh, mikeio.Dfsu2DH))
 
         if ("x" in df) and ("y" in df):
             xy = df[["x", "y"]].to_numpy()
@@ -69,7 +65,7 @@ def validate_positions(cls, mesh, df):
                 "Could not find 'x', 'y' or 'position' columns in DataFrame"
             )
 
-        inside = mesh.contains(xy)
+        inside = mesh.geometry.contains(xy)
         elemid, dist = mesh.geometry.find_nearest_elements(xy, return_distances=True)
         new_positions = mesh.geometry.element_coordinates[elemid, :2]
 

fmdap/pfs_helper.py

@@ -1,14 +1,10 @@
-from collections import namedtuple
+import mikeio
 import pandas as pd
-from mikeio import Pfs
+from mikeio import PfsDocument
 
 
 def pfs2dict(pfs_file):
-    # with open(pfs_file, encoding='iso8859-1') as f:
-    #         pfs = f.read()
-    # y = pfs2yaml(pfs)
-    # d = pfs2yaml.pfs2dict(pfs)
-    d = Pfs(pfs_file).data.to_dict()
+    d = PfsDocument(pfs_file).to_dict()
     return d
 
 
@@ -21,27 +17,9 @@ def get_DA_sections(d):
     return list(get_DA_dict(d).keys())
 
 
-def get_measurements_df(dda):
-    # dda = get_DA_dict(d)
-    meas_sec = dda.get("MEASUREMENTS")
-    if meas_sec is None:
-        raise KeyError("'MEASUREMENTS' section could not be found in dictionary!")
-    n_meas = int(meas_sec.get("number_of_independent_measurements", 0))
+def get_measurements_df(dda: mikeio.PfsSection) -> pd.DataFrame:
+    return dda.MEASUREMENTS.to_dataframe()
 
-    raw = {}
-    for j in range(1, n_meas + 1):
-        raw[j] = meas_sec[f"MEASUREMENT_{j}"]
-    return pd.DataFrame(raw).T
 
-
-def get_diagnostics_df(dda):
-    # dda = get_DA_dict(d)
-    diag_sec = dda.get("DIAGNOSTICS", {}).get("OUTPUTS")
-    if diag_sec is None:
-        raise KeyError("DIAGNOSTICS/OUTPUTS section could not be found in dictionary!")
-    n_diag = int(diag_sec.get("number_of_outputs", 0))
-
-    raw = {}
-    for j in range(1, n_diag + 1):
-        raw[j] = diag_sec[f"OUTPUT_{j}"]
-    return pd.DataFrame(raw).T
+def get_diagnostics_df(dda: mikeio.PfsSection) -> pd.DataFrame:
+    return dda.DIAGNOSTICS.OUTPUTS.to_dataframe()

fmdap/spatial.py

@@ -1,6 +1,6 @@
 import random
 import numpy as np
-from mikeio.spatial.utils import dist_in_meters
+from mikeio.spatial._utils import dist_in_meters
 from scipy.optimize import curve_fit
 
 
@@ -27,12 +27,12 @@ def get_distance_and_corrcoef(dfs, item=0, n_sample=100):
     # assert isinstance(dfs, Dfsu)
 
     elem_ids = None
-    if n_sample is not None and n_sample < dfs.n_elements:
-        n_sample = min(n_sample, dfs.n_elements)
-        elem_ids = random.sample(range(0, dfs.n_elements), n_sample)
+    if n_sample is not None and n_sample < dfs.geometry.n_elements:
+        n_sample = min(n_sample, dfs.geometry.n_elements)
+        elem_ids = random.sample(range(0, dfs.geometry.n_elements), n_sample)
 
-    ec = dfs.element_coordinates[elem_ids, :2]
-    dd = _pairwise_distance(ec, is_geo=dfs.is_geo)
+    ec = dfs.geometry.element_coordinates[elem_ids, :2]
+    dd = _pairwise_distance(ec, is_geo=dfs.geometry.is_geo)
 
     ds = dfs.read(items=[item], elements=elem_ids)[0]
     cc = np.corrcoef(ds.values.T)