Error with "lonlatdepth_dtype=np.float64" in ParticleSet when particles are released with "delay"

[claudiofgcardoso]

Hello,

First of all, I have to thank you for developing this tool and making everything open source! I've been able to learn and fix most issues I encountered through GitHub issues and repositories linked to OceanParcels, and that has been very helpfull :)

I encountered a bug in the lonlatdepth_dtype parameter of ParticleSet.

My Fieldset is a SummedField of ROMS ocean circulation and stokes drift. Sigma depths were already converted to actual depth (m), and horizontal coordinates are on a A-grid. The field is loaded with the following function:

def get_roms(path, beaching = False, start_date=0, finish_date=-1):
    # setting vertical velocities at surface to 0
    def compute(fieldset):
        fieldset.W[0].data[:, 0, :, :] = 0.
        fieldset.W[1].data[:, 0, :, :]=fieldset.W[1].data[:, 0, :, :].clip(min=0)

    fdir=sorted(glob(path + '*.nc'))
    len_path=len(path)   

    if type(start_date) == datetime:
        start_date_str=start_date.strftime("%Y%m%d")
        finish_date_str=finish_date.strftime("%Y%m%d")
        for i, e in enumerate(fdir):
            if start_date != finish_date and e.find(str(start_date_str),len_path) > 0:
                start_index = i 
            elif start_date != finish_date and e.find(str(finish_date_str),len_path) > 0:
                finish_index = i 
            elif start_date == finish_date and e.find(str(start_date_str),len_path) > 0:
                start_index = i 
                finish_index = i 

    if start_date > finish_date:
        nc_filenames = fdir[finish_index:start_index+1]
    else:
        nc_filenames = fdir[start_index:finish_index+1]

    nc_grid=pre_path+'roms_grid_Zlev.nc'

    filenames = {'U': {'lon': nc_grid, 'lat': nc_grid, 'depth': nc_grid, 'data': nc_filenames},
                 'V': {'lon': nc_grid, 'lat': nc_grid, 'depth': nc_grid, 'data': nc_filenames},
                 'W': {'lon': nc_grid, 'lat': nc_grid, 'depth': nc_grid, 'data': nc_filenames}}

    variables_roms = {'U': 'u',
                      'V': 'v',
                      'W': 'w'}

    variables_stokes = {'U': 'u_stokes',
                        'V': 'v_stokes',
                        'W': 'w_stokes'}

    dimensions = {'lon': 'longitude', 
                  'lat': 'latitude', 
                  'depth': 'depth', 
                  'time': 'time'}

    fieldset_roms = FieldSet.from_netcdf(filenames, variables_roms, dimensions, vmax = 10, vmin=-10)
    fieldset_stokes = FieldSet.from_netcdf(filenames, variables_stokes, dimensions, vmax = 10, vmin=-10)

    fieldset_roms.W.set_scaling_factor(-1)
    fieldset_stokes.W.set_scaling_factor(-1)

    fields_W = {'W': fieldset_roms.W + fieldset_stokes.W}
    fields_U = fieldset_roms.U + fieldset_stokes.U
    fields_V = fieldset_roms.V + fieldset_stokes.V

    fieldset_final=FieldSet(fields_U, fields_V, fields_W)    
    fieldset_final.compute_on_defer = compute
    return fieldset_final

My simulation releases a group of 5 particles in different locations from 12 to 12 hours. Following your suggestion of setting the ParticleSet coordinates to double precision in order to prevent uncontrolled particle beaching due to numerical rounding errors, I set lonlatdepth_dtype = np.float64:

def get_particle_set(fieldset, locs):
    class PlasticParticle(JITParticle):
        age = Variable('age', dtype=np.float64, initial=0)
        beached = Variable('beached', dtype=np.float64, initial=0)

    repeatdt = timedelta(hours=12).total_seconds()  # release from the same set of locations every day
    pset = ParticleSet(fieldset, PlasticParticle,
                   lon=locs[:, 1],
                   lat=locs[:, 0],
                   depth=locs[:, 2], 
                   repeatdt = repeatdt, # this is the line where we set a release interval
                   lonlatdepth_dtype=np.float64) # double the coordinate precision to np.float64 to prevent uncontrolled particle beaching, due to numerical rounding errors
    return pset

It did manage to create and save the coordinates of the first set of particles as np.float64, but the second set of particles released 12h later were saved as np.float32. This returns the following error:

INFO: Compiled PlasticParticleAdvectionRK4_beachingBeachTestingAgeing ==> /tmp/parcels-1000/31216d591c72e6060fa8591ebe6953c2.so
 22% (39600.0 of 172800.0) |##           | Elapsed Time: 0:00:21 ETA:   0:01:59Traceback (most recent call last):

  File "<ipython-input-1-52dacf0e2935>", line 429, in <module>
    execute_run(pset, kernel, start_date, finish_date, output_file)

  File "<ipython-input-1-52dacf0e2935>", line 371, in execute_run
    output_file=output_file)

  File "/home/ccardoso/miniconda3/lib/python3.6/site-packages/parcels/particleset.py", line 416, in execute
    self.add(pset_new)

  File "/home/ccardoso/miniconda3/lib/python3.6/site-packages/parcels/particleset.py", line 258, in add
    self._particle_data = np.append(self._particle_data, particles_data)

  File "/home/ccardoso/miniconda3/lib/python3.6/site-packages/numpy/lib/function_base.py", line 4694, in append
    return concatenate((arr, values), axis=axis)

TypeError: invalid type promotion

I believe this happens because the linear interpolation method in the Fieldset (fieldset.U.interp_method = ‘linear') sets the coordinates to np.float32 by default. In order to fix this, I added lonlatdepth_dtype=self.lonlatdepth_dtype to line 410 in particleset.py, and now it works perfectly!

Cheers, Cláudio

[delandmeterp]

Hi Cláudio, Thanks for your interest, your kind words and of course for pointing out and solving this issue! We'll fix it asap!

Cheers, Philippe

[delandmeterp]

We've just checked the code. This problem has actually been fixed in PR #569. (It is that one you mean, right?) Maybe you were using conda's Parcels v2.0.0b2 and not the last master version? We'll release a newer version 2.0 very soon!

[claudiofgcardoso]

Hi Philippe,

Thank you for the replies! You are right, I'm using the V2.0.0b2 because of the issue reported in #547 .. I see that there is a new master version in conda forge a few hours ago, I will install it!

Cheers, Cláudio