Repository: OceanParcels/parcels
Pull request #208 (2017-07-05 15:27:56)

Same arguments for recovery kernels as for normal kernels
This PR adds the same arguments `(particle, fieldset, time, dt)` to RecoveryKernels as to normal Kernels, making it easier to explain how kernels should be implemented

This PR also add pretty-printing of time when `time_origin != 0` using the datetime library


Patch:
diff --git a/examples/example_globcurrent.py b/examples/example_globcurrent.py
index eddf705f..a5f611c5 100644
--- a/examples/example_globcurrent.py
+++ b/examples/example_globcurrent.py
@@ -78,3 +78,15 @@ def test_globcurrent_particles(mode):
 
     assert(abs(pset[0].lon - 23.8) < 1)
     assert(abs(pset[0].lat - -35.3) < 1)
+
+
+@pytest.mark.xfail(reason="Time extrapolation error expected to be thrown")
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_globcurrent_time_extrapolation_error(mode):
+    fieldset = set_globcurrent_fieldset()
+
+    pset = ParticleSet(fieldset, pclass=ptype[mode], lon=[25], lat=[-35])
+    starttime = fieldset.U.time[0]-delta(days=1).total_seconds()
+
+    pset.execute(AdvectionRK4, starttime=starttime, runtime=delta(days=1),
+                 dt=delta(minutes=5), interval=delta(hours=1))
diff --git a/examples/example_recursive_errorhandling.py b/examples/example_recursive_errorhandling.py
index ab36a763..af186cc7 100644
--- a/examples/example_recursive_errorhandling.py
+++ b/examples/example_recursive_errorhandling.py
@@ -41,7 +41,7 @@ def test_recursive_errorhandling(mode, xdim=2, ydim=2):
         if particle.lon <= fieldset.minlon:
             return ErrorCode.Error
 
-    def Error_RandomiseLon(particle):
+    def Error_RandomiseLon(particle, fieldset, time, dt):
         """Error handling kernel that draws a new longitude.
         Note that this new longitude can be smaller than fieldset.minlon"""
         particle.lon = random.uniform(0., 1.)
diff --git a/examples/tutorial_Agulhasparticles.ipynb b/examples/tutorial_Agulhasparticles.ipynb
index 42ae6539..cbc6a2d5 100644
--- a/examples/tutorial_Agulhasparticles.ipynb
+++ b/examples/tutorial_Agulhasparticles.ipynb
@@ -2,21 +2,30 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "## Tutorial showing how to create Parcels in Agulhas animated gif"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "This brief tutorial shows how to recreate the [animated gif](http://oceanparcels.org/animated-gifs/globcurrent_fullyseeded.gif) showing particles in the Agulhas region south of Africa."
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "We start with importing the relevant modules"
    ]
@@ -25,7 +34,9 @@
    "cell_type": "code",
    "execution_count": 1,
    "metadata": {
-    "collapsed": true
+    "collapsed": true,
+    "deletable": true,
+    "editable": true
    },
    "outputs": [],
    "source": [
@@ -36,7 +47,10 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "Now load the Globcurrent fields from the `GlobCurrent_example_data` directory"
    ]
@@ -44,7 +58,11 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false,
+    "deletable": true,
+    "editable": true
+   },
    "outputs": [
     {
      "name": "stderr",
@@ -69,7 +87,10 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "Now create vectors of Longitude and Latitude starting locations on a regular mesh, and use these to initialise a `ParticleSet` object."
    ]
@@ -77,7 +98,11 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {},
+   "metadata": {
+    "collapsed": true,
+    "deletable": true,
+    "editable": true
+   },
    "outputs": [],
    "source": [
     "lons, lats = np.meshgrid(range(15, 35), range(-40, -30))\n",
@@ -86,7 +111,10 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "Now we want to advect the particles. However, the Globcurrent data that we loaded in is only for a limited, regional domain and particles might be able to leave this domain. We therefore need to tell Parcels that particles that leave the domain need to be deleted. We do that using a `Recovery Kernel`, which will be invoked when a particle encounters an `ErrorOutOfBounds` error:"
    ]
@@ -95,17 +123,22 @@
    "cell_type": "code",
    "execution_count": 4,
    "metadata": {
-    "collapsed": true
+    "collapsed": true,
+    "deletable": true,
+    "editable": true
    },
    "outputs": [],
    "source": [
-    "def DeleteParticle(particle):\n",
+    "def DeleteParticle(particle, fieldset, time, dt):\n",
     "    particle.delete()"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "Now we can advect the particles. Note that we do this inside a `for`-loop, so we can save a plot every six hours (which is the value of `runtime`). See the [plotting tutorial](http://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/examples/tutorial_plotting.ipynb) for more information on the `pset.show()` method."
    ]
@@ -113,30 +146,16 @@
   {
    "cell_type": "code",
    "execution_count": 5,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false,
+    "deletable": true,
+    "editable": true
+   },
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py:1767: MatplotlibDeprecationWarning: The get_axis_bgcolor function was deprecated in version 2.0. Use get_facecolor instead.\n",
-      "  axisbgc = ax.get_axis_bgcolor()\n",
-      "/Users/erik/Codes/PARCELScode/parcels/particleset.py:402: RuntimeWarning: invalid value encountered in divide\n",
-      "  normU = U/speed\n",
-      "/Users/erik/Codes/PARCELScode/parcels/particleset.py:403: RuntimeWarning: invalid value encountered in divide\n",
-      "  normV = V/speed\n",
-      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py:3707: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.\n",
-      "  b = ax.ishold()\n",
-      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py:3716: MatplotlibDeprecationWarning: axes.hold is deprecated.\n",
-      "    See the API Changes document (http://matplotlib.org/api/api_changes.html)\n",
-      "    for more details.\n",
-      "  ax.hold(b)\n",
-      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py:3222: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.\n",
-      "  b = ax.ishold()\n",
-      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py:3231: MatplotlibDeprecationWarning: axes.hold is deprecated.\n",
-      "    See the API Changes document (http://matplotlib.org/api/api_changes.html)\n",
-      "    for more details.\n",
-      "  ax.hold(b)\n",
       "INFO: Plot saved to particles00.png\n",
       "INFO: Compiled JITParticleAdvectionRK4 ==> /var/folders/r2/8593q8z93kd7t4j9kbb_f7p00000gr/T/parcels-504/27805ff3aa34ba12ddb373f3f2cb1d1b.so\n",
       "INFO: Plot saved to particles01.png\n",
@@ -159,14 +178,20 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "This now has created 3 plots. Note that the original animated gif contained 20 plots, but to keep running of this notebook fast we have reduced the number here. Of course, it is trivial to increase the number of plots by changing the value in the `range()` in the cell above."
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "As a final step, you can use [ImageMagick](http://www.imagemagick.org/script/index.php) or an online tool to stitch these individual plots together in an animated gif."
    ]
diff --git a/parcels/field.py b/parcels/field.py
index e132ef7a..442faa58 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -32,11 +32,10 @@ class TimeExtrapolationError(RuntimeError):
     """Utility error class to propagate erroneous time extrapolation sampling"""
 
     def __init__(self, time, field=None):
-        self.field = field
-        self.time = time
-        message = "%s sampled outside time domain at time %f." % (
-            field.name if field else "Field", self.time
-        )
+        if field is not None and field.time_origin != 0:
+            time = field.time_origin + timedelta(seconds=time)
+        message = "%s sampled outside time domain at time %s." % (
+            field.name if field else "Field", time)
         message += " Try setting allow_time_extrapolation to True"
         super(TimeExtrapolationError, self).__init__(message)
 
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 80018c44..0139d3d5 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -190,7 +190,7 @@ class Kernel(object):
             for p in error_particles:
                 recovery_kernel = recovery_map[p.state]
                 p.state = ErrorCode.Success
-                recovery_kernel(p)
+                recovery_kernel(p, self.fieldset, p.time, dt)
 
             # Remove all particles that signalled deletion
             remove_deleted(pset)
diff --git a/parcels/kernels/error.py b/parcels/kernels/error.py
index 87c0238d..47546bec 100644
--- a/parcels/kernels/error.py
+++ b/parcels/kernels/error.py
@@ -1,5 +1,6 @@
 """Collection of pre-built recovery kernels"""
 from enum import IntEnum
+from datetime import timedelta
 
 
 __all__ = ['ErrorCode', 'KernelError', 'OutOfBoundsError', 'recovery_map']
@@ -17,61 +18,69 @@ class ErrorCode(IntEnum):
 class KernelError(RuntimeError):
     """General particle kernel error with optional custom message"""
 
-    def __init__(self, particle, msg=None):
-        message = ("%s\nParticle %s\nTime time: %f,\ttimestep size: %f\n") % (
-            particle.state, particle, particle.time, particle.dt
+    def __init__(self, particle, fieldset=None, msg=None):
+        message = ("%s\nParticle %s\nTime: %s,\ttimestep dt: %f\n") % (
+            particle.state, particle, parse_particletime(particle.time, fieldset),
+            particle.dt
         )
         if msg:
             message += msg
         super(KernelError, self).__init__(message)
 
 
-def recovery_kernel_error(particle):
+def parse_particletime(time, fieldset):
+    if fieldset is not None and fieldset.U.time_origin != 0:
+        # TODO assuming that error was thrown on U field
+        time = fieldset.U.time_origin + timedelta(seconds=time)
+    return time
+
+
+def recovery_kernel_error(particle, fieldset, time, dt):
     """Default error kernel that throws exception"""
     msg = "Error: %s" % particle.exception if particle.exception else None
-    raise KernelError(particle, msg=msg)
+    raise KernelError(particle, fieldset=fieldset, msg=msg)
 
 
 class OutOfBoundsError(KernelError):
     """Particle kernel error for out-of-bounds field sampling"""
 
-    def __init__(self, particle, lon=None, lat=None, depth=None, field=None):
+    def __init__(self, particle, fieldset=None, lon=None, lat=None, depth=None):
         if lon and lat:
-            message = "%s sampled at (%f, %f, %f)" % (
-                field.name if field else "Field", lon, lat, depth
+            message = "Field sampled at (%f, %f, %f)" % (
+                lon, lat, depth
             )
         else:
             message = "Out-of-bounds sampling by particle at (%f, %f, %f)" % (
                 particle.lon, particle.lat, particle.depth
             )
-        super(OutOfBoundsError, self).__init__(particle, msg=message)
+        super(OutOfBoundsError, self).__init__(particle, fieldset=fieldset, msg=message)
 
 
 class OutOfTimeError(KernelError):
     """Particle kernel error for time extrapolation field sampling"""
 
-    def __init__(self, particle):
-        message = "Field sampled outside time domain at time %f." % (
-            particle.time
+    def __init__(self, particle, fieldset):
+        message = "Field sampled outside time domain at time %s." % (
+            parse_particletime(particle.time, fieldset)
         )
         message += " Try setting allow_time_extrapolation to True"
-        super(OutOfTimeError, self).__init__(particle, msg=message)
+        super(OutOfTimeError, self).__init__(particle, fieldset=fieldset, msg=message)
 
 
-def recovery_kernel_out_of_bounds(particle):
+def recovery_kernel_out_of_bounds(particle, fieldset, time, dt):
     """Default sampling error kernel that throws OutOfBoundsError"""
     if particle.exception is None:
         # TODO: JIT does not yet provide the context that created
         # the exception. We need to pass that info back from C.
-        raise OutOfBoundsError(particle)
+        raise OutOfBoundsError(particle, fieldset)
     else:
         error = particle.exception
-        raise OutOfBoundsError(particle, error.x, error.y, error.z, error.field)
+        raise OutOfBoundsError(particle, fieldset, error.x, error.y, error.z)
 
 
-def recovery_kernel_time_extrapolation(particle):
+def recovery_kernel_time_extrapolation(particle, fieldset, time, dt):
     """Default sampling error kernel that throws OutOfTimeError"""
-    raise OutOfTimeError(particle)
+    raise OutOfTimeError(particle, fieldset)
 
 
 # Default mapping of failure types (KernelOp)


Test patch:
diff --git a/tests/test_kernel_execution.py b/tests/test_kernel_execution.py
index ced4857b..f6447b00 100644
--- a/tests/test_kernel_execution.py
+++ b/tests/test_kernel_execution.py
@@ -128,7 +128,7 @@ def test_execution_recover_out_of_bounds(fieldset, mode, npart=2):
         fieldset.U[time, particle.lon + 0.1, particle.lat, particle.depth]
         particle.lon += 0.1
 
-    def MoveLeft(particle):
+    def MoveLeft(particle, fieldset, time, dt):
         particle.lon -= 1.
 
     lon = np.linspace(0.05, 0.95, npart, dtype=np.float32)
@@ -147,7 +147,7 @@ def test_execution_delete_out_of_bounds(fieldset, mode, npart=10):
         fieldset.U[time, particle.lon + 0.1, particle.lat, particle.depth]
         particle.lon += 0.1
 
-    def DeleteMe(particle):
+    def DeleteMe(particle, fieldset, time, dt):
         particle.delete()
 
     lon = np.linspace(0.05, 0.95, npart, dtype=np.float32)