Merge pull request #208 from OceanParcels/time_extrapolation_datetime_error

Same arguments for recovery kernels as for normal kernels. Note that after this API change, all recovery kernels will need to be changed in code

Author: erikvansebille

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))

examples/example_recursive_errorhandling.py

@@ -41,7 +41,7 @@ def TestLon(particle, fieldset, time, dt):
         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.)

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,15 +47,22 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "deletable": true,
+    "editable": true
+   },
    "source": [
     "Now load the Globcurrent fields from the `GlobCurrent_example_data` directory"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false,
+    "deletable": true,
+    "editable": true
+   },
    "outputs": [
     {
      "name": "stderr",
@@ -69,15 +87,22 @@
   },
   {
    "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."
    ]
   },
   {
    "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,48 +123,39 @@
    "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."
    ]
   },
   {
    "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."
    ]

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)
 

parcels/kernel.py

@@ -190,7 +190,7 @@ def remove_deleted(pset):
             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)

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)

tests/test_kernel_execution.py

@@ -128,7 +128,7 @@ def MoveRight(particle, fieldset, time, dt):
         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 MoveRight(particle, fieldset, time, dt):
         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)