fixed a bug in RKF time step accounting.

The set tolerances were not always respected.

glacier_flow_tools/pathlines.py

@@ -204,14 +204,13 @@ def compute_pathline(
     time = np.empty(0, dtype=float)
     error_estimate = np.empty(0, dtype=float)
 
-    pts = np.vstack([pts, x])
+    # pts = np.vstack([pts, x])
     vel = f(point, start_time, *f_args)
     v = np.sqrt(vel[0] ** 2 + vel[1] ** 2)
-    velocities = np.vstack([velocities, vel])
-    time = np.append(time, start_time)
-    error_estimate = np.append(error_estimate, 0.0)
+    # velocities = np.vstack([velocities, vel])
+    # time = np.append(time, start_time)
+    # error_estimate = np.append(error_estimate, 0.0)
 
-    k = 0
     p_bar = tqdm_notebook if notebook else tqdm_script
     with p_bar(desc="Integrating pathline", total=end_time, **progress_kwargs) if progress else nullcontext():
         while (t < end_time) and (v > v_threshold):
@@ -240,25 +239,24 @@ def compute_pathline(
 
             r = norm(r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6) / h
             if r <= tol:
+
+                pts = np.append(pts, [x], axis=0)
+                velocities = np.append(velocities, [vel], axis=0)
+                time = np.append(time, t)
+                error_estimate = np.append(error_estimate, r)
+
                 t = t + h
                 x = x + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
 
+                vel = f(x, t, *f_args)
+                v = np.sqrt(vel[0] ** 2 + vel[1] ** 2)
+
             s = (tol / r) ** 0.25
             h = np.minimum(h * s, hmax)
 
             if (h < hmin) and (t < end_time):
-                raise RuntimeError(
-                    f"Error: Could not converge to the required tolerance {tol:e} with minimum stepsize  {hmin:e}"
-                )
-
-            vel = f(point, start_time, *f_args)
-            v = np.sqrt(vel[0] ** 2 + vel[1] ** 2)
-
-            pts = np.append(pts, [x], axis=0)
-            velocities = np.append(velocities, [vel], axis=0)
-            time = np.append(time, t)
-            error_estimate = np.append(error_estimate, r)
-            k += 1
+                print(f"Error: Could not converge to the required tolerance {tol:e} with minimum stepsize  {hmin:e}")
+                continue
 
     return pts, velocities, time, error_estimate
 

glacier_flow_tools/profiles.py

@@ -193,9 +193,9 @@ def plot_obs_sims_profile(
 
     if interactive:
         # The standard backend is not thread-safe, but 'agg' works with the dask client.
-        matplotlib.use("agg")
-    else:
         matplotlib.use("module://matplotlib_inline.backend_inline")
+    else:
+        matplotlib.use("agg")
 
     fig = ds.profiles.plot_obs_sims(
         sigma=sigma,
@@ -1309,10 +1309,10 @@ def plot(
         """
 
         if interactive:
+            matplotlib.use("module://matplotlib_inline.backend_inline")
+        else:
             # The standard backend is not thread-safe, but 'agg' works with the dask client.
             matplotlib.use("agg")
-        else:
-            matplotlib.use("module://matplotlib_inline.backend_inline")
 
         plt.rcParams["font.size"] = fontsize
         n_exps = self._obj["exp_id"].size

notebooks/create_starting_points.ipynb

@@ -0,0 +1,102 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4d8a8e39-8429-43db-b71d-62c38611922f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geopandas as gp\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import shapely\n",
+    "from shapely import Point"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "08335445-fffc-4641-bb97-cb2161270636",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "glaciers = gp.read_file(\"../glacier_flow_tools/data/greenland-flux-gates-5.gpkg\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "03a00076-8cbc-48dd-9670-2fb831d0fcb9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "distance = 20_000 # m\n",
+    "spacing = 4_000 # m"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3400ba1c-1a35-4182-9279-17784006143e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%time\n",
+    "glaciers_starting_points = []\n",
+    "for _, glacier in glaciers.iterrows():\n",
+    "    x_c, y_c = glacier.geometry.centroid.coords[0]\n",
+    "    x = np.linspace(x_c - distance, x_c + distance, int(distance * 2 / spacing + 1))\n",
+    "    y = np.linspace(y_c - distance, y_c + distance, int(distance * 2 / spacing + 1))\n",
+    "    X, Y = np.meshgrid(x, y)\n",
+    "\n",
+    "    gl = [glacier.copy() for _ in range(len(X.ravel()))]\n",
+    "    for m, (x, y) in enumerate(zip(X.ravel(), Y.ravel())):\n",
+    "        gl[m][\"pt\"] = m\n",
+    "        gl[m].geometry = Point(x, y)\n",
+    "\n",
+    "    glaciers_starting_points.append(gp.GeoDataFrame(gl))\n",
+    "\n",
+    "glaciers_starting_points = pd.concat(glaciers_starting_points)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c5f9024b-c6d8-4437-bc06-24872ce57a18",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "glaciers_starting_points.to_file(f\"starting_pts-5-{spacing}m.gpkg\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f781801-5e72-46f9-95e9-675ab947c2dc",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

pathlines/compute_pathlines.py

@@ -20,6 +20,7 @@
 Calculate pathlines (trajectories).
 """
 
+import time
 from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
 from pathlib import Path
 
@@ -58,7 +59,7 @@
         type=float,
         default=1.0,
     )
-    parser.add_argument("--tol", help="""Adaptive time stepping tolerance.""", type=float, default=1.0)
+    parser.add_argument("--tol", help="""Adaptive time stepping tolerance. Default=1e-3""", type=float, default=1e-3)
     parser.add_argument(
         "--start_time",
         help="""Start time. Default=0.0""",
@@ -77,6 +78,13 @@
         action="store_true",
         default=False,
     )
+    parser.add_argument(
+        "--output_type",
+        help="""Save result as Points or LineStrings.""",
+        choices=["point", "line"],
+        default="point",
+    )
+
     parser.add_argument(
         "--v_threshold",
         help="""Threshold velocity below which solver stops Default is 0.0.""",
@@ -89,8 +97,6 @@
 
     p = Path(options.outfile[-1])
     p.parent.mkdir(parents=True, exist_ok=True)
-    line_p = Path("line_" + options.outfile[-1])
-    line_p.parent.mkdir(parents=True, exist_ok=True)
 
     df = gp.read_file(options.vector_url)
     starting_points_df = geopandas_dataframe_shorten_lines(df).convert.to_points()
@@ -108,6 +114,7 @@
 
     n_pts = len(starting_points_df)
 
+    start = time.time()
     with tqdm_joblib(
         tqdm(desc="Processing Pathlines", total=n_pts, leave=True, position=0)
     ) as progress_bar:  # pylint: disable=unused-variable
@@ -123,21 +130,23 @@
                 end_time=options.end_time,
                 v_threshold=options.v_threshold,
                 progress=False,
-                progress_kwargs={"leave": False, "position": index},
             )
             for index, df in starting_points_df.iterrows()
         )
-    ps = [
-        series_to_pathline_geopandas_dataframe(s.drop("geometry", errors="ignore"), pathlines[k])
-        for k, s in starting_points_df.iterrows()
-    ]
+    time_elapsed = time.time() - start
+    print(f"Time elapsed {time_elapsed:.0f}s\n")
 
-    result = pd.concat(ps).reset_index(drop=True)
-    result.to_file(p, mode="w")
+    print(f"Saving {p}")
 
-    ps = [
-        pathline_to_line_geopandas_dataframe(pathlines[k][0], attrs={"pathline_id": [k]})
-        for k, _ in starting_points_df.iterrows()
-    ]
+    if options.output_type == "point":
+        ps = [
+            series_to_pathline_geopandas_dataframe(s.drop("geometry", errors="ignore"), pathlines[k])
+            for k, s in starting_points_df.iterrows()
+        ]
+    else:
+        ps = [
+            pathline_to_line_geopandas_dataframe(pathlines[k][0], attrs={"pathline_id": [k]})
+            for k, _ in starting_points_df.iterrows()
+        ]
     result = pd.concat(ps).reset_index(drop=True)
-    result.to_file(line_p, mode="w")
+    result.to_file(p, mode="w")

profiles/compute_profiles.py

@@ -326,5 +326,5 @@ def __call__(self, parser, namespace, values, option_string=None):
             progress(futures_computed)
             futures_gathered = client.gather(futures_computed)
 
-    time_elapsed = time.time() - start
-    print(f"Time elapsed {time_elapsed:.0f}s")
+            time_elapsed = time.time() - start
+            print(f"Time elapsed {time_elapsed:.0f}s")