Repository: OceanParcels/parcels
Pull request #569 (2019-03-20 08:02:20)

Double precision particlefile netcdf
Adding support for writing of float64 variables in netcdf files


Patch:
diff --git a/parcels/particlefile.py b/parcels/particlefile.py
index da5b5c9f..61775be2 100644
--- a/parcels/particlefile.py
+++ b/parcels/particlefile.py
@@ -87,19 +87,24 @@ class ParticleFile(object):
             self.time.calendar = 'standard' if self.particleset.time_origin.calendar == 'np_datetime64' else self.particleset.time_origin.calendar
         self.time.axis = "T"
 
-        self.lat = self.dataset.createVariable("lat", "f4", coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        if self.particleset.lonlatdepth_dtype is np.float64:
+            lonlatdepth_precision = "f8"
+        else:
+            lonlatdepth_precision = "f4"
+
+        self.lat = self.dataset.createVariable("lat", lonlatdepth_precision, coords, fill_value=np.nan, chunksizes=self.chunksizes)
         self.lat.long_name = ""
         self.lat.standard_name = "latitude"
         self.lat.units = "degrees_north"
         self.lat.axis = "Y"
 
-        self.lon = self.dataset.createVariable("lon", "f4", coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        self.lon = self.dataset.createVariable("lon", lonlatdepth_precision, coords, fill_value=np.nan, chunksizes=self.chunksizes)
         self.lon.long_name = ""
         self.lon.standard_name = "longitude"
         self.lon.units = "degrees_east"
         self.lon.axis = "X"
 
-        self.z = self.dataset.createVariable("z", "f4", coords, fill_value=np.nan, chunksizes=self.chunksizes)
+        self.z = self.dataset.createVariable("z", lonlatdepth_precision, coords, fill_value=np.nan, chunksizes=self.chunksizes)
         self.z.long_name = ""
         self.z.standard_name = "depth"
         self.z.units = "m"
diff --git a/parcels/particleset.py b/parcels/particleset.py
index 8ae4564e..a8973414 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -413,7 +413,8 @@ class ParticleSet(object):
             if abs(time-next_prelease) < tol:
                 pset_new = ParticleSet(fieldset=self.fieldset, time=time, lon=self.repeatlon,
                                        lat=self.repeatlat, depth=self.repeatdepth,
-                                       pclass=self.repeatpclass, **self.repeatkwargs)
+                                       pclass=self.repeatpclass, lonlatdepth_dtype=self.lonlatdepth_dtype,
+                                       **self.repeatkwargs)
                 for p in pset_new:
                     p.dt = dt
                 self.add(pset_new)


Test patch:
diff --git a/tests/test_particle_sets.py b/tests/test_particle_sets.py
index ef2d7b88..160f9c91 100644
--- a/tests/test_particle_sets.py
+++ b/tests/test_particle_sets.py
@@ -417,6 +417,21 @@ def test_variable_written_once(fieldset, mode, tmpdir, npart):
     assert [v == 0 for v in vfile]
 
 
+@pytest.mark.parametrize('mode', ['scipy', 'jit'])
+def test_variable_write_double(fieldset, mode, tmpdir):
+    filepath = tmpdir.join("pfile_variable_write_double")
+
+    def Update_lon(particle, fieldset, time):
+        particle.lon += 0.1
+
+    pset = ParticleSet(fieldset, pclass=JITParticle, lon=[0], lat=[0], lonlatdepth_dtype=np.float64)
+    pset.execute(pset.Kernel(Update_lon), endtime=1, dt=0.1,
+                 output_file=pset.ParticleFile(name=filepath, outputdt=0.1))
+    ncfile = Dataset(filepath+".nc", 'r', 'NETCDF4')
+    lons = ncfile.variables['lon'][:]
+    assert (isinstance(lons[0, 0], np.float64))
+
+
 @pytest.mark.parametrize('mode', ['scipy', 'jit'])
 def test_variable_written_ondelete(fieldset, mode, tmpdir, npart=3):
     filepath = tmpdir.join("pfile_on_delete_written_variables")