Repository: OceanParcels/parcels
Pull request #496 (2018-11-26 08:12:43)

Pset from field


Patch:
diff --git a/parcels/particleset.py b/parcels/particleset.py
index fa305ed4..55eb8704 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -150,7 +150,7 @@ class ParticleSet(object):
         :param fieldset: :mod:`parcels.fieldset.FieldSet` object from which to sample velocity
         :param pclass: mod:`parcels.particle.JITParticle` or :mod:`parcels.particle.ScipyParticle`
                  object that defines custom particle
-        :param start_field: Field for initialising particles stochastically according to the presented density field.
+        :param start_field: Field for initialising particles stochastically (horizontally)  according to the presented density field.
         :param size: Initial size of particle set
         :param mode: Type of random sampling. Currently only 'monte_carlo' is implemented
         :param depth: Optional list of initial depth values for particles. Default is 0m
@@ -159,7 +159,15 @@ class ParticleSet(object):
         """
 
         if mode == 'monte_carlo':
-            p_interior = np.squeeze(start_field.data[:, :-1, :-1])
+            if start_field.interp_method == 'cgrid_tracer':
+                p_interior = np.squeeze(start_field.data[0, 1:, 1:])
+            else:  # if A-grid
+                d = start_field.data
+                p_interior = (d[0, :-1, :-1] + d[0, 1:, :-1] + d[0, :-1, 1:] + d[0, 1:, 1:])/4.
+                p_interior = np.where(d[0, :-1, :-1] == 0, 0, p_interior)
+                p_interior = np.where(d[0, 1:, :-1] == 0, 0, p_interior)
+                p_interior = np.where(d[0, 1:, 1:] == 0, 0, p_interior)
+                p_interior = np.where(d[0, :-1, 1:] == 0, 0, p_interior)
             p = np.reshape(p_interior, (1, p_interior.size))
             inds = np.random.choice(p_interior.size, size, replace=True, p=p[0] / np.sum(p))
             xsi = np.random.uniform(size=len(inds))
diff --git a/parcels/plotting.py b/parcels/plotting.py
index b9362d0f..be467849 100644
--- a/parcels/plotting.py
+++ b/parcels/plotting.py
@@ -159,13 +159,14 @@ def plotfield(field, show_time=None, domain=None, depth_level=0, projection=None
         if field[0].interp_method == 'cgrid_velocity':
             logger.warning_once('Plotting a C-grid velocity field is achieved via an A-grid projection, reducing the plot accuracy')
             d = np.empty_like(data[0])
-            d[:-1, :] = (data[:-1, :] + data[1:, :]) / 2.
-            d[-1, :] = data[-1, :]
+            d[:-1, :] = (data[0][:-1, :] + data[0][1:, :]) / 2.
+            d[-1, :] = data[0][-1, :]
             data[0] = d
             d = np.empty_like(data[0])
-            d[:, :-1] = (data[:, :-1] + data[:, 1:]) / 2.
-            d[:, -1] = data[:, -1]
+            d[:, :-1] = (data[0][:, :-1] + data[0][:, 1:]) / 2.
+            d[:, -1] = data[0][:, -1]
             data[1] = d
+
         spd = data[0] ** 2 + data[1] ** 2
         speed = np.sqrt(spd, where=spd > 0)
         vmin = speed.min() if vmin is None else vmin


Test patch:
diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py
index 43dd5a93..dbbfa3a5 100644
--- a/tests/test_particle_sets.py
+++ b/tests/test_particle_sets.py
@@ -430,3 +430,52 @@ def test_variable_written_ondelete(fieldset, mode, tmpdir, npart=3):
     ncfile = Dataset(filepath+".nc", 'r', 'NETCDF4')
     lon = ncfile.variables['lon'][:]
     assert (lon.size == noutside)
+
+
+@pytest.mark.parametrize('staggered_grid', ['Agrid', 'Cgrid'])
+def test_from_field_exact_val(staggered_grid):
+    xdim = 4
+    ydim = 3
+
+    lon = np.linspace(-1, 2, xdim, dtype=np.float32)
+    lat = np.linspace(50, 52, ydim, dtype=np.float32)
+
+    dimensions = {'lat': lat, 'lon': lon}
+    if staggered_grid == 'Agrid':
+        U = np.zeros((ydim, xdim), dtype=np.float32)
+        V = np.zeros((ydim, xdim), dtype=np.float32)
+        data = {'U': np.array(U, dtype=np.float32), 'V': np.array(V, dtype=np.float32)}
+        mask = np.array([[1, 1, 0, 0],
+                         [1, 1, 1, 0],
+                         [1, 1, 1, 1]])
+        fieldset = FieldSet.from_data(data, dimensions, mesh='flat')
+
+        FMask = Field('mask', mask, lon, lat)
+        fieldset.add_field(FMask)
+    elif staggered_grid == 'Cgrid':
+        U = np.array([[0, 0, 0, 0],
+                      [1, 0, 0, 0],
+                      [1, 1, 0, 0]])
+        V = np.array([[0, 1, 0, 0],
+                      [0, 1, 0, 0],
+                      [0, 1, 1, 0]])
+        data = {'U': np.array(U, dtype=np.float32), 'V': np.array(V, dtype=np.float32)}
+        mask = np.array([[-1, -1, -1, -1],
+                         [-1, 1, 0, 0],
+                         [-1, 1, 1, 0]])
+        fieldset = FieldSet.from_data(data, dimensions, mesh='flat')
+        fieldset.U.interp_method = 'cgrid_velocity'
+        fieldset.V.interp_method = 'cgrid_velocity'
+
+        FMask = Field('mask', mask, lon, lat, interp_method='cgrid_tracer')
+        fieldset.add_field(FMask)
+
+    class SampleParticle(ptype['scipy']):
+        mask = Variable('mask', initial=fieldset.mask)
+
+    pset = ParticleSet.from_field(fieldset, size=400, pclass=SampleParticle, start_field=FMask, time=0)
+    # pset.show(field=FMask)
+    assert np.allclose([p.mask for p in pset], 1)
+    assert (np.array([p.lon for p in pset]) <= 1).all()
+    test = np.logical_or(np.array([p.lon for p in pset]) <= 0, np.array([p.lat for p in pset]) >= 51)
+    assert test.all()