Repository: OceanParcels/parcels
Pull request #310 (2018-03-08 22:25:59)

Reducing file size of ParticleFile when particles are deleted and not all started at same time
This PR reduces the file set of ParticleFile, by making the time index where each particle Variable is written dependent on the Particle. This means that if Particles are added during `.execute()`, they are still written from index 0, reducing the file size if particles are also deleted.

Since for each Particle, we also store `time`, we still store all information required for post-processing.

For example, in the function `test_pset_repeated_release_delayed_adding_deleting(fieldset(), 'jit', repeatdt=2, tmpdir='', dt=1, maxvar=4, runtime=10)`, we go from
```
ncfile.variables['time'][:] =
[[0.0 1.0 2.0 3.0 4.0 5.0 -- -- -- -- --]
 [-- -- 2.0 3.0 4.0 5.0 6.0 7.0 -- -- --]
 [-- -- -- -- 4.0 5.0 6.0 7.0 8.0 9.0 --]
 [-- -- -- -- -- -- 6.0 7.0 8.0 9.0 10.0]
 [-- -- -- -- -- -- -- -- 8.0 9.0 10.0]
 [-- -- -- -- -- -- -- -- -- -- 10.0]]
```
To
```
ncfile.variables['time'][:] =
[[0.0 1.0 2.0 3.0 4.0 5.0]
 [2.0 3.0 4.0 5.0 6.0 7.0]
 [4.0 5.0 6.0 7.0 8.0 9.0]
 [6.0 7.0 8.0 9.0 10.0 --]
 [8.0 9.0 10.0 -- -- --]
 [10.0 -- -- -- -- --]]
```



Patch:
diff --git a/parcels/particlefile.py b/parcels/particlefile.py
index 3c887922..24cb1e49 100644
--- a/parcels/particlefile.py
+++ b/parcels/particlefile.py
@@ -102,7 +102,7 @@ class ParticleFile(object):
                 getattr(self, v.name).standard_name = v.name
                 getattr(self, v.name).units = "unknown"
 
-        self.idx = 0
+        self.idx = np.empty(shape=0)
 
     def __del__(self):
         self.dataset.close()
@@ -125,25 +125,25 @@ class ParticleFile(object):
            (self.write_ondelete is False or deleted_only is True):
             if pset.size > 0:
 
-                first_write = [p for p in pset if p.fileid < 0]
+                first_write = [p for p in pset if p.fileid < 0 or len(self.idx) == 0]  # len(self.idx)==0 in case pset is written to new ParticleFile
                 for p in first_write:
                     p.fileid = self.lasttraj
                     self.lasttraj += 1
 
-                inds = [p.fileid for p in pset]
-
-                self.id[inds, self.idx] = [p.id for p in pset]
-                self.time[inds, self.idx] = time
-                self.lat[inds, self.idx] = np.array([p.lat for p in pset])
-                self.lon[inds, self.idx] = np.array([p.lon for p in pset])
-                self.z[inds, self.idx] = np.array([p.depth for p in pset])
-                for var in self.user_vars:
-                    getattr(self, var)[inds, self.idx] = np.array([getattr(p, var) for p in pset])
-                for var in self.user_vars_once:
-                    if np.any(first_write):
-                        vals = [getattr(p, var) for p in first_write]
-                        newinds = [p.fileid for p in first_write]
-                        getattr(self, var)[newinds] = np.array(vals)
+                self.idx = np.append(self.idx, np.zeros(len(first_write)))
+
+                for p in pset:
+                    i = p.fileid
+                    self.id[i, self.idx[i]] = p.id
+                    self.time[i, self.idx[i]] = time
+                    self.lat[i, self.idx[i]] = p.lat
+                    self.lon[i, self.idx[i]] = p.lon
+                    self.z[i, self.idx[i]] = p.depth
+                    for var in self.user_vars:
+                        getattr(self, var)[i, self.idx[i]] = getattr(p, var)
+                for p in first_write:
+                    for var in self.user_vars_once:
+                        getattr(self, var)[p.fileid] = getattr(p, var)
             else:
                 logger.warning("ParticleSet is empty on writing as array")
 
diff --git a/parcels/scripts/plotParticles.py b/parcels/scripts/plotParticles.py
index 9338f07f..9917c171 100644
--- a/parcels/scripts/plotParticles.py
+++ b/parcels/scripts/plotParticles.py
@@ -32,9 +32,10 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
         return
 
     pfile = Dataset(filename, 'r')
-    lon = pfile.variables['lon']
-    lat = pfile.variables['lat']
-    z = pfile.variables['z']
+    lon = np.ma.filled(pfile.variables['lon'], np.nan)
+    lat = np.ma.filled(pfile.variables['lat'], np.nan)
+    time = np.ma.filled(pfile.variables['time'], np.nan)
+    z = np.ma.filled(pfile.variables['z'], np.nan)
 
     if(recordedvar is not None):
         record = pfile.variables[recordedvar]
@@ -62,13 +63,19 @@ def plotTrajectoriesFile(filename, mode='2d', tracerfile=None, tracerfield='P',
     elif mode == 'movie2d' or 'movie2d_notebook':
         fig = plt.figure()
         ax = plt.axes(xlim=(np.nanmin(lon), np.nanmax(lon)), ylim=(np.nanmin(lat), np.nanmax(lat)))
-        scat = ax.scatter(lon[:, 0], lat[:, 0], s=60, cmap=plt.get_cmap('autumn'))  # cmaps not working?
-        frames = np.arange(1, lon.shape[1])
+        plottimes = np.unique(time)
+        plottimes = plottimes[~np.isnan(plottimes)]
+        b = time == plottimes[0]
+        scat = ax.scatter(lon[b], lat[b], s=60, cmap=plt.get_cmap('autumn'))  # cmaps not working?
+        ttl = ax.set_title('Particle at time ' + str(plottimes[0]))
+        frames = np.arange(1, len(plottimes))
 
         def animate(t):
-            scat.set_offsets(np.matrix((lon[:, t], lat[:, t])).transpose())
+            b = time == plottimes[t]
+            scat.set_offsets(np.matrix((lon[b], lat[b])).transpose())
+            ttl.set_text('Particle at time ' + str(plottimes[t]))
             if recordedvar is not None:
-                scat.set_array(record[:, t])
+                scat.set_array(record[b])
             return scat,
 
         rc('animation', html='html5')


Test patch:
diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py
index 8dbf92aa..d03fcffc 100644
--- a/tests/test_particle_sets.py
+++ b/tests/test_particle_sets.py
@@ -101,7 +101,7 @@ def test_pset_repeated_release_delayed_adding_deleting(fieldset, mode, repeatdt,
     assert np.allclose([p.sample_var for p in pset], np.arange(maxvar, -1, -repeatdt))
     ncfile = Dataset(outfilepath+".nc", 'r', 'NETCDF4')
     samplevar = ncfile.variables['sample_var'][:]
-    assert samplevar.shape == (runtime // repeatdt+1, runtime+1)
+    assert samplevar.shape == (runtime // repeatdt+1, min(maxvar+1, runtime)+1)
     if repeatdt == 0:
         # test whether samplevar[i, i+k] = k for k=range(maxvar)
         for k in range(maxvar):