Merge branch 'master' into summedvectorfield_eval_implementation

Author: erikvansebille

parcels/compilation/codecompiler.py

@@ -237,7 +237,10 @@ def _create_compile_process_(self, cmd, src, log):
 Source/Destination file: %s
 Log file: %s
 
-Log output: %s""" % (" ".join(cmd), src, logfile.name, logfile2.read())
+Log output: %s
+If you are on macOS, it might help to type 'export CC=gcc'
+
+""" % (" ".join(cmd), src, logfile.name, logfile2.read())
                 raise RuntimeError(err)
         return True
 

parcels/compilation/codegenerator.py

@@ -331,6 +331,10 @@ def visit_Attribute(self, node):
                 raise NotImplementedError("Cannot convert numpy functions in kernels to C-code.\n"
                                           "Either use functions from the math library or run Parcels in Scipy mode.\n"
                                           "For more information, see http://oceanparcels.org/faq.html#kernelwriting")
+            elif node.value.id in ['random']:
+                raise NotImplementedError("Cannot convert random functions in kernels to C-code.\n"
+                                          "Use `import parcels.rng as ParcelsRandom` and then ParcelsRandom.random(), ParcelsRandom.uniform() etc.\n"
+                                          "For more information, see http://oceanparcels.org/faq.html#kernelwriting")
             else:
                 raise NotImplementedError("Cannot convert '%s' used in kernel to C-code" % node.value.id)
 
@@ -692,9 +696,9 @@ def visit_If(self, node):
 
     def visit_Compare(self, node):
         self.visit(node.left)
-        assert(len(node.ops) == 1)
+        assert (len(node.ops) == 1)
         self.visit(node.ops[0])
-        assert(len(node.comparators) == 1)
+        assert (len(node.comparators) == 1)
         self.visit(node.comparators[0])
         node.ccode = "%s %s %s" % (node.left.ccode, node.ops[0].ccode,
                                    node.comparators[0].ccode)

parcels/examples/documentation_indexing.ipynb

@@ -65,7 +65,11 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {},
+   "metadata": {
+    "tags": [
+     "raises-exception"
+    ]
+   },
    "outputs": [
     {
      "ename": "ValueError",
@@ -123,19 +127,7 @@
    "cell_type": "code",
    "execution_count": 4,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "WARNING: File NemoNorthSeaORCA025-N006_data/coordinates.nc could not be decoded properly by xarray (version 0.15.1).\n",
-      "         It will be opened with no decoding. Filling values might be wrongly parsed.\n",
-      "WARNING: Casting lon data to np.float32\n",
-      "WARNING: Casting lat data to np.float32\n",
-      "WARNING: Trying to initialize a shared grid with different chunking sizes - action prohibited. Replacing requested field_chunksize with grid's master chunksize.\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "mesh_mask = data_path + 'coordinates.nc'\n",
     "\n",
@@ -147,7 +139,14 @@
     "              'V': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'},\n",
     "              'W': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'}}\n",
     "\n",
-    "fieldsetC = FieldSet.from_nemo(filenames, variables, dimensions)"
+    "fieldsetC = FieldSet.from_nemo(filenames, variables, dimensions, netcdf_decodewarning=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note by the way, that we used `netcdf_decodewarning=False` in the `FieldSet.from_nemo()` call above. This is to silence an expected warning because the time dimension in the `coordinates.nc` file can't be decoded by `xarray`."
    ]
   },
   {
@@ -221,7 +220,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -235,7 +234,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.8.13"
   }
  },
  "nbformat": 4,

parcels/examples/documentation_unstuck_Agrid.ipynb

@@ -11,7 +11,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In [another notebook](https://ndViewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_stuck_particles.ipynb), we have shown how particles may end up getting stuck on land, especially in A gridded velocity fields. Here we show how you can work around this problem and how large the effects of the solutions on the trajectories are.\n",
+    "In [another notebook](https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_stuck_particles.ipynb), we have shown how particles may end up getting stuck on land, especially in A gridded velocity fields. Here we show how you can work around this problem and how large the effects of the solutions on the trajectories are.\n",
     "\n",
     "Common solutions are:\n",
     "1. [Delete the particles](#1.-Particle-deletion)\n",
@@ -28,7 +28,7 @@
     "1. Flag particles within a specific distance to the shore\n",
     "2. Flag particles in any gridcell that has a shore edge\n",
     "\n",
-    "As argued in the [previous notebook](https://ndViewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_stuck_particles.ipynb), it is important to accurately plot the grid discretization, in order to understand the motion of particles near the boundary. The velocity fields can best be depicted using points or arrows that define the velocity at a single position. Four of these nodes then form gridcells that can be shown using tiles, for example with `matplotlib.pyplot.pcolormesh`."
+    "As argued in the [previous notebook](https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_stuck_particles.ipynb), it is important to accurately plot the grid discretization, in order to understand the motion of particles near the boundary. The velocity fields can best be depicted using points or arrows that define the velocity at a single position. Four of these nodes then form gridcells that can be shown using tiles, for example with `matplotlib.pyplot.pcolormesh`."
    ]
   },
   {

parcels/examples/example_dask_chunk_OCMs.py

@@ -281,8 +281,8 @@ def test_globcurrent_2D(mode, chunk_mode):
     elif chunk_mode == 'failsafe':  # chunking time but not lat
         assert (len(field_set.U.grid.load_chunk) != 1)
         assert (len(field_set.V.grid.load_chunk) != 1)
-    assert(abs(pset[0].lon - 23.8) < 1)
-    assert(abs(pset[0].lat - -35.3) < 1)
+    assert abs(pset[0].lon - 23.8) < 1
+    assert abs(pset[0].lat - -35.3) < 1
 
 
 @pytest.mark.parametrize('mode', ['jit'])
@@ -382,8 +382,8 @@ def test_ofam_3D(mode, chunk_mode):
         matching_fields = (field_set.U.grid.chunk_info == field_set.V.grid.chunk_info)
         matching_uniformblocks = (len(field_set.U.grid.load_chunk) == (1 * int(math.ceil(1.0/60.0)) * int(math.ceil(601.0/50.0)) * int(math.ceil(2001.0/100.0))))
         assert (matching_uniformblocks or (matching_fields and matching_numblocks))
-    assert(abs(pset[0].lon - 173) < 1)
-    assert(abs(pset[0].lat - 11) < 1)
+    assert abs(pset[0].lon - 173) < 1
+    assert abs(pset[0].lat - 11) < 1
 
 
 @pytest.mark.parametrize('mode', ['jit'])

parcels/examples/example_globcurrent.py

@@ -39,10 +39,10 @@ def set_globcurrent_fieldset(filename=None, indices=None, deferred_load=True, us
 @pytest.mark.parametrize('use_xarray', [True, False])
 def test_globcurrent_fieldset(use_xarray):
     fieldset = set_globcurrent_fieldset(use_xarray=use_xarray)
-    assert(fieldset.U.lon.size == 81)
-    assert(fieldset.U.lat.size == 41)
-    assert(fieldset.V.lon.size == 81)
-    assert(fieldset.V.lat.size == 41)
+    assert fieldset.U.lon.size == 81
+    assert fieldset.U.lat.size == 41
+    assert fieldset.V.lon.size == 81
+    assert fieldset.V.lat.size == 41
 
     if not use_xarray:
         indices = {'lon': [5], 'lat': range(20, 30)}
@@ -88,8 +88,8 @@ def test_globcurrent_particles(mode, use_xarray):
 
     pset.execute(AdvectionRK4, runtime=delta(days=1), dt=delta(minutes=5))
 
-    assert(abs(pset[0].lon - 23.8) < 1)
-    assert(abs(pset[0].lat - -35.3) < 1)
+    assert abs(pset[0].lon - 23.8) < 1
+    assert abs(pset[0].lat - -35.3) < 1
 
 
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])

parcels/examples/example_moving_eddies.py

@@ -157,8 +157,8 @@ def test_moving_eddies_fieldset(mode, mesh, tmpdir):
         assert (pset[0].lon < 2.2e5 and 1.1e5 < pset[0].lat < 1.2e5)
         assert (pset[1].lon < 2.2e5 and 3.7e5 < pset[1].lat < 3.8e5)
     else:
-        assert(pset[0].lon < 2.0 and 46.2 < pset[0].lat < 46.25)
-        assert(pset[1].lon < 2.0 and 48.8 < pset[1].lat < 48.85)
+        assert (pset[0].lon < 2.0 and 46.2 < pset[0].lat < 46.25)
+        assert (pset[1].lon < 2.0 and 48.8 < pset[1].lat < 48.85)
 
 
 def fieldsetfile(mesh, tmpdir):

parcels/examples/example_ofam.py

@@ -33,7 +33,7 @@ def set_ofam_fieldset(deferred_load=True, use_xarray=False):
 def test_ofam_fieldset_fillvalues(use_xarray):
     fieldset = set_ofam_fieldset(deferred_load=False, use_xarray=use_xarray)
     # V.data[0, 0, 150] is a landpoint, that makes NetCDF4 generate a masked array, instead of an ndarray
-    assert(fieldset.V.data[0, 0, 150] == 0)
+    assert fieldset.V.data[0, 0, 150] == 0
 
 
 @pytest.mark.parametrize('dt', [delta(minutes=-5), delta(minutes=5)])
@@ -66,5 +66,5 @@ def test_ofam_particles(mode, use_xarray):
 
     pset.execute(AdvectionRK4, runtime=delta(days=10), dt=delta(minutes=5))
 
-    assert(abs(pset[0].lon - 173) < 1)
-    assert(abs(pset[0].lat - 11) < 1)
+    assert abs(pset[0].lon - 173) < 1
+    assert abs(pset[0].lat - 11) < 1

parcels/examples/example_peninsula.py

@@ -145,10 +145,10 @@ def test_peninsula_fieldset(mode, mesh, tmpdir):
     pset = peninsula_example(fieldset, outfile, 5, mode=mode, degree=1)
     # Test advection accuracy by comparing streamline values
     err_adv = np.abs(pset.p_start - pset.p)
-    assert(err_adv <= 1.e-3).all()
+    assert (err_adv <= 1.e-3).all()
     # Test Field sampling accuracy by comparing kernel against Field sampling
     err_smpl = np.array([abs(pset.p[i] - pset.fieldset.P[0., pset.depth[i], pset.lat[i], pset.lon[i]]) for i in range(pset.size)])
-    assert(err_smpl <= 1.e-3).all()
+    assert (err_smpl <= 1.e-3).all()
 
 
 @pytest.mark.parametrize('mode', ['scipy'])  # Analytical Advection only implemented in Scipy mode
@@ -161,7 +161,7 @@ def test_peninsula_fieldset_AnalyticalAdvection(mode, mesh, tmpdir):
                              method=AdvectionAnalytical)
     # Test advection accuracy by comparing streamline values
     err_adv = np.array([abs(p.p_start - p.p) for p in pset])
-    assert(err_adv <= 1.e-1).all()
+    assert (err_adv <= 1.e-1).all()
 
 
 def fieldsetfile(mesh, tmpdir):
@@ -182,10 +182,10 @@ def test_peninsula_file(mode, mesh, tmpdir):
     pset = peninsula_example(fieldset, outfile, 5, mode=mode, degree=1)
     # Test advection accuracy by comparing streamline values
     err_adv = np.abs(pset.p_start - pset.p)
-    assert(err_adv <= 1.e-3).all()
+    assert (err_adv <= 1.e-3).all()
     # Test Field sampling accuracy by comparing kernel against Field sampling
     err_smpl = np.array([abs(pset.p[i] - pset.fieldset.P[0., pset.depth[i], pset.lat[i], pset.lon[i]]) for i in range(pset.size)])
-    assert(err_smpl <= 1.e-3).all()
+    assert (err_smpl <= 1.e-3).all()
 
 
 if __name__ == "__main__":

parcels/examples/example_radial_rotation.py

@@ -31,11 +31,11 @@ def radial_rotation_fieldset(xdim=200, ydim=200):  # Define 2D flat, square fiel
         for j in range(lat.size):
 
             r = np.sqrt((lon[i]-x0)**2 + (lat[j]-y0)**2)  # Define radial displacement.
-            assert(r >= 0.)
-            assert(r <= np.sqrt(x0**2 + y0**2))
+            assert r >= 0.
+            assert r <= np.sqrt(x0**2 + y0**2)
 
             theta = math.atan2((lat[j]-y0), (lon[i]-x0))  # Define the polar angle.
-            assert(abs(theta) <= np.pi)
+            assert abs(theta) <= np.pi
 
             U[j, i] = r * math.sin(theta) * omega
             V[j, i] = -r * math.cos(theta) * omega
@@ -75,10 +75,10 @@ def test_rotation_example(mode, tmpdir):
     fieldset = radial_rotation_fieldset()
     outfile = tmpdir.join("RadialParticle")
     pset = rotation_example(fieldset, outfile, mode=mode)
-    assert(pset[0].lon == 30. and pset[0].lat == 30.)  # Particle at centre of Field remains stationary.
+    assert (pset[0].lon == 30. and pset[0].lat == 30.)  # Particle at centre of Field remains stationary.
     vals = true_values(pset.time[1])
-    assert(np.allclose(pset[1].lon, vals[0], 1e-5))    # Check advected values against calculated values.
-    assert(np.allclose(pset[1].lat, vals[1], 1e-5))
+    assert np.allclose(pset[1].lon, vals[0], 1e-5)    # Check advected values against calculated values.
+    assert np.allclose(pset[1].lat, vals[1], 1e-5)
 
 
 if __name__ == "__main__":