Time float, increase of amplitude precision

CHANGELOG.rst

@@ -10,6 +10,11 @@ and this project adheres to `Semantic Versioning <https://semver.org/spec/v2.0.0
 ------------
 Changed
 ^^^^^^^
+- Now time will be allow second precision in storage on uint32 from 01/01/1950 to 01/01/2086
+  new identification will be produce with this type, old file could be still loaded.
+  If you use old identification to track use `--unraw` option to unpack old time and store in new format.
+- Now amplitude is stored with .1 mm of precision, same advice than time.
+
 Fixed
 ^^^^^
 - GridCollection get_next_time_step & get_previous_time_step needed more files to work in the dataset list.

doc/run_tracking.rst

@@ -26,7 +26,7 @@ Example of conf.yaml
 
     # Number of timestep for missing detection
     VIRTUAL_LENGTH_MAX: 3
-    # Minimal time to consider as a full track
+    # Minimal number of timesteps to consider as a long track
     TRACK_DURATION_MIN: 10
 
 To run:
@@ -69,7 +69,7 @@ With yaml you could also select another tracker:
 
     # Number of timesteps for missing detection
     VIRTUAL_LENGTH_MAX: 3
-    # Minimal time to consider as a full track
+    # Minimal number of timesteps to consider as a long track
     TRACK_DURATION_MIN: 10
 
     CLASS:

examples/06_grid_manipulation/pet_okubo_weiss.py

@@ -1,6 +1,6 @@
 r"""
 Get Okubo Weis
-=====================
+==============
 
 .. math:: OW = S_n^2 + S_s^2 + \omega^2
 

examples/16_network/pet_follow_particle.py

@@ -133,12 +133,12 @@ def update(frame):
 shape = (n.obs.size, 2)
 # Forward run
 i_target_f, pct_target_f = -ones(shape, dtype="i4"), zeros(shape, dtype="i1")
-for t in range(t_start, t_end - dt):
+for t in arange(t_start, t_end - dt):
     particle_candidate(c, n, step, t, i_target_f, pct_target_f, n_days=dt)
 
 # Backward run
 i_target_b, pct_target_b = -ones(shape, dtype="i4"), zeros(shape, dtype="i1")
-for t in range(t_start + dt, t_end):
+for t in arange(t_start + dt, t_end):
     particle_candidate(c, n, step, t, i_target_b, pct_target_b, n_days=-dt)
 
 # %%
@@ -149,6 +149,11 @@ def update(frame):
 ax_2nd_f = fig.add_axes([0.52, 0.05, 0.45, 0.45])
 ax_1st_b.set_title("Backward advection for each time step")
 ax_1st_f.set_title("Forward advection for each time step")
+ax_1st_b.set_ylabel("Color -> First target\nLatitude")
+ax_2nd_b.set_ylabel("Color -> Secondary target\nLatitude")
+ax_2nd_b.set_xlabel("Julian days"), ax_2nd_f.set_xlabel("Julian days")
+ax_1st_f.set_yticks([]), ax_2nd_f.set_yticks([])
+ax_1st_f.set_xticks([]), ax_1st_b.set_xticks([])
 
 
 def color_alpha(target, pct, vmin=5, vmax=80):

notebooks/python_module/16_network/pet_follow_particle.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\n# Follow particle\n"
+        "\nFollow particle\n===============\n"
       ]
     },
     {
@@ -55,7 +55,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Schema\n\n"
+        "Schema\n------\n\n"
       ]
     },
     {
@@ -73,7 +73,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Animation\nParticle settings\n\n"
+        "Animation\n---------\nParticle settings\n\n"
       ]
     },
     {
@@ -109,7 +109,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "### Particle advection\n\n"
+        "Particle advection\n^^^^^^^^^^^^^^^^^^\n\n"
       ]
     },
     {
@@ -120,7 +120,7 @@
       },
       "outputs": [],
       "source": [
-        "step = 1 / 60.0\n\nx, y = meshgrid(arange(24, 36, step), arange(31, 36, step))\nx0, y0 = x.reshape(-1), y.reshape(-1)\n# Pre-order to speed up\n_, i = group_obs(x0, y0, 1, 360)\nx0, y0 = x0[i], y0[i]\n\nt_start, t_end = n.period\ndt = 14\n\nshape = (n.obs.size, 2)\n# Forward run\ni_target_f, pct_target_f = -ones(shape, dtype=\"i4\"), zeros(shape, dtype=\"i1\")\nfor t in range(t_start, t_end - dt):\n    particle_candidate(x0, y0, c, n, t, i_target_f, pct_target_f, n_days=dt)\n\n# Backward run\ni_target_b, pct_target_b = -ones(shape, dtype=\"i4\"), zeros(shape, dtype=\"i1\")\nfor t in range(t_start + dt, t_end):\n    particle_candidate(x0, y0, c, n, t, i_target_b, pct_target_b, n_days=-dt)"
+        "step = 1 / 60.0\n\nx, y = meshgrid(arange(24, 36, step), arange(31, 36, step))\nx0, y0 = x.reshape(-1), y.reshape(-1)\n# Pre-order to speed up\n_, i = group_obs(x0, y0, 1, 360)\nx0, y0 = x0[i], y0[i]\n\nt_start, t_end = n.period\ndt = 14\n\nshape = (n.obs.size, 2)\n# Forward run\ni_target_f, pct_target_f = -ones(shape, dtype=\"i4\"), zeros(shape, dtype=\"i1\")\nfor t in arange(t_start, t_end - dt):\n    particle_candidate(x0, y0, c, n, t, i_target_f, pct_target_f, n_days=dt)\n\n# Backward run\ni_target_b, pct_target_b = -ones(shape, dtype=\"i4\"), zeros(shape, dtype=\"i1\")\nfor t in arange(t_start + dt, t_end):\n    particle_candidate(x0, y0, c, n, t, i_target_b, pct_target_b, n_days=-dt)"
       ]
     },
     {
@@ -131,7 +131,7 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure(figsize=(10, 10))\nax_1st_b = fig.add_axes([0.05, 0.52, 0.45, 0.45])\nax_2nd_b = fig.add_axes([0.05, 0.05, 0.45, 0.45])\nax_1st_f = fig.add_axes([0.52, 0.52, 0.45, 0.45])\nax_2nd_f = fig.add_axes([0.52, 0.05, 0.45, 0.45])\nax_1st_b.set_title(\"Backward advection for each time step\")\nax_1st_f.set_title(\"Forward advection for each time step\")\n\n\ndef color_alpha(target, pct, vmin=5, vmax=80):\n    color = cmap(n.segment[target])\n    # We will hide under 5 % and from 80% to 100 % it will be 1\n    alpha = (pct - vmin) / (vmax - vmin)\n    alpha[alpha < 0] = 0\n    alpha[alpha > 1] = 1\n    color[:, 3] = alpha\n    return color\n\n\nkw = dict(\n    name=None, yfield=\"longitude\", event=False, zorder=-100, s=(n.speed_area / 20e6)\n)\nn.scatter_timeline(ax_1st_b, c=color_alpha(i_target_b.T[0], pct_target_b.T[0]), **kw)\nn.scatter_timeline(ax_2nd_b, c=color_alpha(i_target_b.T[1], pct_target_b.T[1]), **kw)\nn.scatter_timeline(ax_1st_f, c=color_alpha(i_target_f.T[0], pct_target_f.T[0]), **kw)\nn.scatter_timeline(ax_2nd_f, c=color_alpha(i_target_f.T[1], pct_target_f.T[1]), **kw)\nfor ax in (ax_1st_b, ax_2nd_b, ax_1st_f, ax_2nd_f):\n    n.display_timeline(ax, field=\"longitude\", marker=\"+\", lw=2, markersize=5)\n    ax.grid()"
+        "fig = plt.figure(figsize=(10, 10))\nax_1st_b = fig.add_axes([0.05, 0.52, 0.45, 0.45])\nax_2nd_b = fig.add_axes([0.05, 0.05, 0.45, 0.45])\nax_1st_f = fig.add_axes([0.52, 0.52, 0.45, 0.45])\nax_2nd_f = fig.add_axes([0.52, 0.05, 0.45, 0.45])\nax_1st_b.set_title(\"Backward advection for each time step\")\nax_1st_f.set_title(\"Forward advection for each time step\")\nax_1st_b.set_ylabel(\"Color -> First target\\nLatitude\")\nax_2nd_b.set_ylabel(\"Color -> Secondary target\\nLatitude\")\nax_2nd_b.set_xlabel(\"Julian days\"), ax_2nd_f.set_xlabel(\"Julian days\")\nax_1st_f.set_yticks([]), ax_2nd_f.set_yticks([])\nax_1st_f.set_xticks([]), ax_1st_b.set_xticks([])\n\n\ndef color_alpha(target, pct, vmin=5, vmax=80):\n    color = cmap(n.segment[target])\n    # We will hide under 5 % and from 80% to 100 % it will be 1\n    alpha = (pct - vmin) / (vmax - vmin)\n    alpha[alpha < 0] = 0\n    alpha[alpha > 1] = 1\n    color[:, 3] = alpha\n    return color\n\n\nkw = dict(\n    name=None, yfield=\"longitude\", event=False, zorder=-100, s=(n.speed_area / 20e6)\n)\nn.scatter_timeline(ax_1st_b, c=color_alpha(i_target_b.T[0], pct_target_b.T[0]), **kw)\nn.scatter_timeline(ax_2nd_b, c=color_alpha(i_target_b.T[1], pct_target_b.T[1]), **kw)\nn.scatter_timeline(ax_1st_f, c=color_alpha(i_target_f.T[0], pct_target_f.T[0]), **kw)\nn.scatter_timeline(ax_2nd_f, c=color_alpha(i_target_f.T[1], pct_target_f.T[1]), **kw)\nfor ax in (ax_1st_b, ax_2nd_b, ax_1st_f, ax_2nd_f):\n    n.display_timeline(ax, field=\"longitude\", marker=\"+\", lw=2, markersize=5)\n    ax.grid()"
       ]
     }
   ],
@@ -151,7 +151,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.9"
+      "version": "3.7.7"
     }
   },
   "nbformat": 4,

notebooks/python_module/16_network/pet_segmentation_anim.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\n# Network segmentation process\n"
+        "\nNetwork segmentation process\n============================\n"
       ]
     },
     {
@@ -62,7 +62,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Load data\nLoad data where observations are put in same network but no segmentation\n\n"
+        "Load data\n---------\nLoad data where observations are put in same network but no segmentation\n\n"
       ]
     },
     {
@@ -80,7 +80,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Do segmentation\nSegmentation based on maximum overlap, temporal window for candidates = 5 days\n\n"
+        "Do segmentation\n---------------\nSegmentation based on maximum overlap, temporal window for candidates = 5 days\n\n"
       ]
     },
     {
@@ -98,7 +98,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Anim\n\n"
+        "Anim\n----\n\n"
       ]
     },
     {
@@ -109,14 +109,14 @@
       },
       "outputs": [],
       "source": [
-        "def update(i_frame):\n    tr = TRACKS[i_frame]\n    mappable_tracks.set_array(tr)\n    s = 40 * ones(tr.shape)\n    s[tr == 0] = 4\n    mappable_tracks.set_sizes(s)\n\n    indices_frames = INDICES[i_frame]\n    mappable_CONTOUR.set_data(\n        e.contour_lon_e[indices_frames],\n        e.contour_lat_e[indices_frames],\n    )\n    mappable_CONTOUR.set_color(cmap.colors[tr[indices_frames] % len(cmap.colors)])\n    return (mappable_tracks,)\n\n\nfig = plt.figure(figsize=(16, 9), dpi=60)\nax = fig.add_axes([0.04, 0.06, 0.94, 0.88], projection=GUI_AXES)\nax.set_title(f\"{len(e)} observations to segment\")\nax.set_xlim(19, 29), ax.set_ylim(31, 35.5), ax.grid()\nvmax = TRACKS[-1].max()\ncmap = ListedColormap([\"gray\", *e.COLORS[:-1]], name=\"from_list\", N=vmax)\nmappable_tracks = ax.scatter(\n    e.lon, e.lat, c=TRACKS[0], cmap=cmap, vmin=0, vmax=vmax, s=20\n)\nmappable_CONTOUR = ax.plot(\n    e.contour_lon_e[INDICES[0]], e.contour_lat_e[INDICES[0]], color=cmap.colors[0]\n)[0]\nani = VideoAnimation(fig, update, frames=range(1, len(TRACKS), 4), interval=125)"
+        "def update(i_frame):\n    tr = TRACKS[i_frame]\n    mappable_tracks.set_array(tr)\n    s = 40 * ones(tr.shape)\n    s[tr == 0] = 4\n    mappable_tracks.set_sizes(s)\n\n    indices_frames = INDICES[i_frame]\n    mappable_CONTOUR.set_data(\n        e.contour_lon_e[indices_frames], e.contour_lat_e[indices_frames],\n    )\n    mappable_CONTOUR.set_color(cmap.colors[tr[indices_frames] % len(cmap.colors)])\n    return (mappable_tracks,)\n\n\nfig = plt.figure(figsize=(16, 9), dpi=60)\nax = fig.add_axes([0.04, 0.06, 0.94, 0.88], projection=GUI_AXES)\nax.set_title(f\"{len(e)} observations to segment\")\nax.set_xlim(19, 29), ax.set_ylim(31, 35.5), ax.grid()\nvmax = TRACKS[-1].max()\ncmap = ListedColormap([\"gray\", *e.COLORS[:-1]], name=\"from_list\", N=vmax)\nmappable_tracks = ax.scatter(\n    e.lon, e.lat, c=TRACKS[0], cmap=cmap, vmin=0, vmax=vmax, s=20\n)\nmappable_CONTOUR = ax.plot(\n    e.contour_lon_e[INDICES[0]], e.contour_lat_e[INDICES[0]], color=cmap.colors[0]\n)[0]\nani = VideoAnimation(fig, update, frames=range(1, len(TRACKS), 4), interval=125)"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Final Result\n\n"
+        "Final Result\n------------\n\n"
       ]
     },
     {
@@ -147,7 +147,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.9"
+      "version": "3.7.7"
     }
   },
   "nbformat": 4,

share/tracking.yaml

@@ -4,8 +4,9 @@ PATHS:
   # Path for saving of outputs
   SAVE_DIR: '/home/emason/toto/'
 
-# Minimum number of observations to store eddy
+# Minimal number of timesteps to consider as a long track
 TRACK_DURATION_MIN: 4
+# Number of timesteps for missing detection
 VIRTUAL_LENGTH_MAX: 0
 
 CLASS:

src/py_eddy_tracker/__init__.py

@@ -106,12 +106,17 @@ def parse_args(self, *args, **kwargs):
         return opts
 
 
+TIME_MODELS = ["%Y%m%d", "%Y%m%d%H%M%S", "%Y%m%dT%H%M%S"]
+
+
 VAR_DESCR = dict(
     time=dict(
         attr_name="time",
         nc_name="time",
         old_nc_name=["j1"],
-        nc_type="int32",
+        nc_type="float64",
+        output_type="uint32",
+        scale_factor=1 / 86400.0,
         nc_dims=("obs",),
         nc_attr=dict(
             standard_name="time",
@@ -251,7 +256,7 @@ def parse_args(self, *args, **kwargs):
         old_nc_name=["A"],
         nc_type="float32",
         output_type="uint16",
-        scale_factor=0.001,
+        scale_factor=0.0001,
         nc_dims=("obs",),
         nc_attr=dict(
             long_name="Amplitude",

src/py_eddy_tracker/appli/eddies.py

@@ -15,7 +15,7 @@
 from numpy import bincount, bytes_, empty, in1d, unique
 from yaml import safe_load
 
-from .. import EddyParser
+from .. import TIME_MODELS, EddyParser
 from ..observations.observation import EddiesObservations, reverse_index
 from ..observations.tracking import TrackEddiesObservations
 from ..tracking import Correspondances
@@ -223,7 +223,7 @@ def browse_dataset_in(
     data_dir,
     files_model,
     date_regexp,
-    date_model,
+    date_model=None,
     start_date=None,
     end_date=None,
     sub_sampling_step=1,
@@ -238,11 +238,7 @@ def browse_dataset_in(
         filenames = bytes_(glob(full_path))
 
     dataset_list = empty(
-        len(filenames),
-        dtype=[
-            ("filename", "S500"),
-            ("date", "datetime64[D]"),
-        ],
+        len(filenames), dtype=[("filename", "S500"), ("date", "datetime64[s]")],
     )
     dataset_list["filename"] = filenames
 
@@ -268,10 +264,21 @@ def browse_dataset_in(
                 str_date = result.groups()[0]
 
         if str_date is not None:
-            item["date"] = datetime.strptime(str_date, date_model).date()
+            if date_model is None:
+                model_found = False
+                for model in TIME_MODELS:
+                    try:
+                        item["date"] = datetime.strptime(str_date, model)
+                        model_found = True
+                        break
+                    except ValueError:
+                        pass
+                if not model_found:
+                    raise Exception("No time model found")
+            else:
+                item["date"] = datetime.strptime(str_date, date_model)
 
     dataset_list.sort(order=["date", "filename"])
-
     steps = unique(dataset_list["date"][1:] - dataset_list["date"][:-1])
     if len(steps) > 1:
         raise Exception("Several days steps in grid dataset %s" % steps)
@@ -304,7 +311,7 @@ def track(
     correspondances_only=False,
     **kw_c,
 ):
-    kw = dict(date_regexp=".*_([0-9]*?).[nz].*", date_model="%Y%m%d")
+    kw = dict(date_regexp=".*_([0-9]*?).[nz].*")
     if isinstance(pattern, list):
         kw.update(dict(data_dir=None, files_model=None, files=pattern))
     else:
@@ -323,10 +330,9 @@ def track(
     c = Correspondances(datasets=datasets["filename"], **kw_c)
     c.track()
     logger.info("Track finish")
-    t0, t1 = c.period
     kw_save = dict(
-        date_start=t0,
-        date_stop=t1,
+        date_start=datasets["date"][0],
+        date_stop=datasets["date"][-1],
         date_prod=datetime.now(),
         path=output_dir,
         sign_type=c.current_obs.sign_legend,
@@ -351,11 +357,13 @@ def track(
 
     short_c = c._copy()
     short_c.shorter_than(size_max=nb_obs_min)
-    c.longer_than(size_min=nb_obs_min)
-
-    long_track = c.merge(raw_data=raw)
     short_track = short_c.merge(raw_data=raw)
 
+    if c.longer_than(size_min=nb_obs_min) is False:
+        long_track = short_track.empty_dataset()
+    else:
+        long_track = c.merge(raw_data=raw)
+
     # We flag obs
     if c.virtual:
         long_track["virtual"][:] = long_track["time"] == 0

src/py_eddy_tracker/appli/grid.py

@@ -5,7 +5,7 @@
 from argparse import Action
 from datetime import datetime
 
-from .. import EddyParser
+from .. import TIME_MODELS, EddyParser
 from ..dataset.grid import RegularGridDataset, UnRegularGridDataset
 
 
@@ -121,7 +121,16 @@ def eddy_id(args=None):
         cut_wavelength = [0, *cut_wavelength]
     inf_bnds, upper_bnds = cut_wavelength
 
-    date = datetime.strptime(args.datetime, "%Y%m%d")
+    model_found = False
+    for model in TIME_MODELS:
+        try:
+            date = datetime.strptime(args.datetime, model)
+            model_found = True
+            break
+        except ValueError:
+            pass
+    if not model_found:
+        raise Exception("No time model found")
     kwargs = dict(
         step=args.isoline_step,
         shape_error=args.fit_errmax,
@@ -150,7 +159,7 @@ def eddy_id(args=None):
         sampling_method=args.sampling_method,
         **kwargs,
     )
-    out_name = date.strftime("%(path)s/%(sign_type)s_%Y%m%d.nc")
+    out_name = date.strftime("%(path)s/%(sign_type)s_%Y%m%dT%H%M%S.nc")
     a.write_file(path=args.path_out, filename=out_name, zarr_flag=args.zarr)
     c.write_file(path=args.path_out, filename=out_name, zarr_flag=args.zarr)
 

src/py_eddy_tracker/gui.py

@@ -291,8 +291,8 @@ def get_infos(self, name, index):
         i_first = d.index_from_track[tr]
         track = d.obs[i_first : i_first + nb]
         nb -= 1
-        t0 = timedelta(days=int(track[0]["time"])) + datetime(1950, 1, 1)
-        t1 = timedelta(days=int(track[-1]["time"])) + datetime(1950, 1, 1)
+        t0 = timedelta(days=track[0]["time"]) + datetime(1950, 1, 1)
+        t1 = timedelta(days=track[-1]["time"]) + datetime(1950, 1, 1)
         txt = f"--{name}--\n"
         txt += f"    {t0} -> {t1}\n"
         txt += f"    Tracks : {tr}  {now['n']}/{nb} ({now['n'] / nb * 100:.2f} %)\n"

src/py_eddy_tracker/observations/groups.py

@@ -186,15 +186,15 @@ def filled_by_interpolation(self, mask):
 
         .. minigallery:: py_eddy_tracker.TrackEddiesObservations.filled_by_interpolation
         """
-
+        if self.track.size == 0:
+            return
         nb_filled = mask.sum()
         logger.info("%d obs will be filled (unobserved)", nb_filled)
 
         nb_obs = len(self)
         index = arange(nb_obs)
 
         for field in self.obs.dtype.descr:
-            # print(f"field : {field}")
             var = field[0]
             if (
                 var in ["n", "virtual", "track", "cost_association"]

src/py_eddy_tracker/observations/observation.py

@@ -1328,6 +1328,8 @@ def solve_conflict(cost):
     def solve_simultaneous(cost):
         """Write something (TODO)"""
         mask = ~cost.mask
+        if mask.size == 0:
+            return mask
         # Count number of links by self obs and other obs
         self_links, other_links = sum_row_column(mask)
         max_links = max(self_links.max(), other_links.max())