Merge pull request AntSimi#95 from AntSimi/time_float

Time float, increase of amplitude precision

Author: ludwigVonKoopa

CHANGELOG.rst

@@ -10,13 +10,25 @@ and this project adheres to `Semantic Versioning <https://semver.org/spec/v2.0.0
 ------------
 Changed
 ^^^^^^^
+
+- Now time allows second precision (instead of daily precision) in storage on uint32 from 01/01/1950 to 01/01/2086
+  New identifications are produced with this type, old files could still be loaded.
+  If you use old identifications for tracking use the `--unraw` option to unpack old times and store data with the new format.
+- Now amplitude is stored with .1 mm of precision (instead of 1 mm), same advice as for time.
+
 Fixed
 ^^^^^
-- GridCollection get_next_time_step & get_previous_time_step needed more files to work in the dataset list.
-  The loop needed explicitly self.dataset[i+-1] even when i==0, therefore indice went out of range
+
 Added
 ^^^^^
 
+[3.5.0] - 2021-06-22
+--------------------
+
+Fixed
+^^^^^
+- GridCollection get_next_time_step & get_previous_time_step needed more files to work in the dataset list.
+  The loop needed explicitly self.dataset[i+-1] even when i==0, therefore indice went out of range
 
 [3.4.0] - 2021-03-29
 --------------------

doc/run_tracking.rst

@@ -5,9 +5,9 @@ Tracking
 Requirements
 ************
 
-Before to run tracking, you will need to run identification on every time step of the period (period of your study).
+Before tracking, you will need to run identification on every time step of the period (period of your study).
 
-**Advice** : Before to run tracking, displaying some identification file allows to learn a lot
+**Advice** : Before tracking, displaying some identification files. You will learn a lot
 
 Default method
 **************
@@ -24,9 +24,9 @@ Example of conf.yaml
       FILES_PATTERN: MY_IDENTIFICATION_PATH/Anticyclonic*.nc
       SAVE_DIR: MY_OUTPUT_PATH
 
-    # Number of timestep for missing detection
+    # Number of consecutive timesteps with missing detection allowed
     VIRTUAL_LENGTH_MAX: 3
-    # Minimal time to consider as a full track
+    # Minimal number of timesteps to considered as a long trajectory
     TRACK_DURATION_MIN: 10
 
 To run:
@@ -63,13 +63,13 @@ With yaml you could also select another tracker:
 .. code-block:: yaml
 
     PATHS:
-      # Files produces with EddyIdentification
+      # Files produced with EddyIdentification
       FILES_PATTERN: MY/IDENTIFICATION_PATH/Anticyclonic*.nc
       SAVE_DIR: MY_OUTPUT_PATH
 
-    # Number of timesteps for missing detection
+    # Number of consecutive timesteps with missing detection allowed
     VIRTUAL_LENGTH_MAX: 3
-    # Minimal time to consider as a full track
+    # Minimal number of timesteps to considered as a long trajectory
     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

@@ -1,12 +1,13 @@
 PATHS:
-  # Files produces with EddyIdentification
+  # Files produced with EddyIdentification
   FILES_PATTERN: /home/emason/toto/Anticyclonic_*.nc
-  # Path for saving of outputs
+  # Path to save outputs
   SAVE_DIR: '/home/emason/toto/'
 
-# Minimum number of observations to store eddy
-TRACK_DURATION_MIN: 4
+# Number of consecutive timesteps with missing detection allowed
 VIRTUAL_LENGTH_MAX: 0
+# Minimal number of timesteps to considered as a long trajectory
+TRACK_DURATION_MIN: 4
 
 CLASS:
    MODULE: py_eddy_tracker.featured_tracking.area_tracker

src/py_eddy_tracker/__init__.py

@@ -22,6 +22,7 @@
 
 import logging
 from argparse import ArgumentParser
+from datetime import datetime
 
 import zarr
 
@@ -106,12 +107,26 @@ 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"]
+
+
+def identify_time(str_date):
+    for model in TIME_MODELS:
+        try:
+            return datetime.strptime(str_date, model)
+        except ValueError:
+            pass
+    raise Exception("No time model found")
+
+
 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 +266,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 EddyParser, identify_time
 from ..observations.observation import EddiesObservations, reverse_index
 from ..observations.tracking import TrackEddiesObservations
 from ..tracking import Correspondances
@@ -163,7 +163,12 @@ def eddies_tracking():
     parser.add_argument(
         "--zarr", action="store_true", help="Output will be wrote in zarr"
     )
-    parser.add_argument("--unraw", action="store_true", help="Load unraw data")
+    parser.add_argument(
+        "--unraw",
+        action="store_true",
+        help="Load unraw data, use only for netcdf."
+        "If unraw is active, netcdf is loaded without apply scalefactor and add_offset.",
+    )
     parser.add_argument(
         "--blank_period",
         type=int,
@@ -223,7 +228,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 +243,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,13 +269,15 @@ 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:
+                item["date"] = identify_time(str_date)
+            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)
+        raise Exception("Several timesteps in grid dataset %s" % steps)
 
     if sub_sampling_step != 1:
         logger.info("Grid subsampling %d", sub_sampling_step)
@@ -304,7 +307,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 +326,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 +353,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