Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
216 changes: 1 addition & 215 deletions parcels/particleset.py
Original file line number Diff line number Diff line change
Expand Up @@ -347,221 +347,7 @@ def populate_indices(self):
self._data["ei"][:, i] = idx # assumes that we are in the surface layer (zi=0)

@classmethod
def from_list(
cls, fieldset, pclass, lon, lat, depth=None, time=None, repeatdt=None, lonlatdepth_dtype=None, **kwargs
):
"""Initialise the ParticleSet from lists of lon and lat.

Parameters
----------
fieldset :
mod:`parcels.fieldset.FieldSet` object from which to sample velocity
pclass :
Particle class. May be a parcels.particle.Particle class as defined in parcels, or a subclass defining a custom particle.
lon :
List of initial longitude values for particles
lat :
List of initial latitude values for particles
depth :
Optional list of initial depth values for particles. Default is 0m
time :
Optional list of start time values for particles. Default is fieldset.U.time[0]
repeatdt :
Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
lonlatdepth_dtype :
Floating precision for lon, lat, depth particle coordinates.
It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
and np.float64 if the interpolation method is 'cgrid_velocity'
Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
**kwargs :
Keyword arguments passed to the particleset constructor.
"""
return cls(
fieldset=fieldset,
pclass=pclass,
lon=lon,
lat=lat,
depth=depth,
time=time,
repeatdt=repeatdt,
lonlatdepth_dtype=lonlatdepth_dtype,
**kwargs,
)

@classmethod
def from_line(
cls,
fieldset,
pclass,
start,
finish,
size,
depth=None,
time=None,
repeatdt=None,
lonlatdepth_dtype=None,
**kwargs,
):
"""Create a particleset in the shape of a line (according to a cartesian grid).

Initialise the ParticleSet from start/finish coordinates with equidistant spacing
Note that this method uses simple numpy.linspace calls and does not take into account
great circles, so may not be a exact on a globe

Parameters
----------
fieldset :
mod:`parcels.fieldset.FieldSet` object from which to sample velocity
pclass :
Particle class. May be a parcels.particle.Particle as defined in parcels, or a subclass defining a custom particle.
start :
Start point (longitude, latitude) for initialisation of particles on a straight line.
finish :
End point (longitude, latitude) for initialisation of particles on a straight line.
size :
Initial size of particle set
depth :
Optional list of initial depth values for particles. Default is 0m
time :
Optional start time value for particles. Default is fieldset.U.time[0]
repeatdt :
Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
lonlatdepth_dtype :
Floating precision for lon, lat, depth particle coordinates.
It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
and np.float64 if the interpolation method is 'cgrid_velocity'
"""
lon = np.linspace(start[0], finish[0], size)
lat = np.linspace(start[1], finish[1], size)
if type(depth) in [int, float]:
depth = [depth] * size
return cls(
fieldset=fieldset,
pclass=pclass,
lon=lon,
lat=lat,
depth=depth,
time=time,
repeatdt=repeatdt,
lonlatdepth_dtype=lonlatdepth_dtype,
**kwargs,
)

@classmethod
def _monte_carlo_sample(cls, start_field, size, mode="monte_carlo"):
"""Converts a starting field into a monte-carlo sample of lons and lats.

Parameters
----------
start_field : parcels.field.Field
mod:`parcels.fieldset.Field` object for initialising particles stochastically (horizontally) according to the presented density field.
size :

mode :
(Default value = 'monte_carlo')

Returns
-------
list of float
A list of longitude values.
list of float
A list of latitude values.
"""
if mode == "monte_carlo":
data = start_field.data if isinstance(start_field.data, np.ndarray) else np.array(start_field.data)
if start_field.interp_method == "cgrid_tracer":
p_interior = np.squeeze(data[0, 1:, 1:])
else: # if A-grid
d = data
p_interior = (d[0, :-1, :-1] + d[0, 1:, :-1] + d[0, :-1, 1:] + d[0, 1:, 1:]) / 4.0
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))
eta = np.random.uniform(size=len(inds))
j, i = np.unravel_index(inds, p_interior.shape)
grid = start_field.grid
lon, lat = ([], [])
if grid._gtype in [GridType.RectilinearZGrid, GridType.RectilinearSGrid]:
lon = grid.lon[i] + xsi * (grid.lon[i + 1] - grid.lon[i])
lat = grid.lat[j] + eta * (grid.lat[j + 1] - grid.lat[j])
else:
lons = np.array([grid.lon[j, i], grid.lon[j, i + 1], grid.lon[j + 1, i + 1], grid.lon[j + 1, i]])
if grid.mesh == "spherical":
lons[1:] = np.where(lons[1:] - lons[0] > 180, lons[1:] - 360, lons[1:])
lons[1:] = np.where(-lons[1:] + lons[0] > 180, lons[1:] + 360, lons[1:])
lon = (
(1 - xsi) * (1 - eta) * lons[0]
+ xsi * (1 - eta) * lons[1]
+ xsi * eta * lons[2]
+ (1 - xsi) * eta * lons[3]
)
lat = (
(1 - xsi) * (1 - eta) * grid.lat[j, i]
+ xsi * (1 - eta) * grid.lat[j, i + 1]
+ xsi * eta * grid.lat[j + 1, i + 1]
+ (1 - xsi) * eta * grid.lat[j + 1, i]
)
return list(lat), list(lon)
else:
raise NotImplementedError(f'Mode {mode} not implemented. Please use "monte carlo" algorithm instead.')

@classmethod
def from_field(
cls,
fieldset,
pclass,
start_field,
size,
mode="monte_carlo",
depth=None,
time=None,
repeatdt=None,
lonlatdepth_dtype=None,
):
"""Initialise the ParticleSet randomly drawn according to distribution from a field.

Parameters
----------
fieldset : parcels.fieldset.FieldSet
mod:`parcels.fieldset.FieldSet` object from which to sample velocity
pclass :
Particle class. May be a parcels.particle.Particle class as defined in parcels, or a subclass defining a custom particle.
start_field : parcels.field.Field
Field for initialising particles stochastically (horizontally) according to the presented density field.
size :
Initial size of particle set
mode :
Type of random sampling. Currently only 'monte_carlo' is implemented (Default value = 'monte_carlo')
depth :
Optional list of initial depth values for particles. Default is 0m
time :
Optional start time value for particles. Default is fieldset.U.time[0]
repeatdt :
Optional interval (in seconds) on which to repeat the release of the ParticleSet (Default value = None)
lonlatdepth_dtype :
Floating precision for lon, lat, depth particle coordinates.
It is either np.float32 or np.float64. Default is np.float32 if fieldset.U.interp_method is 'linear'
and np.float64 if the interpolation method is 'cgrid_velocity'
"""
lat, lon = cls._monte_carlo_sample(start_field, size, mode)

return cls(
fieldset=fieldset,
pclass=pclass,
lon=lon,
lat=lat,
depth=depth,
time=time,
lonlatdepth_dtype=lonlatdepth_dtype,
repeatdt=repeatdt,
)

@classmethod
def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime=None, repeatdt=None, **kwargs):
def from_particlefile(cls, fieldset, pclass, filename, restart=True, restarttime=None, **kwargs):
"""Initialise the ParticleSet from a zarr ParticleFile.
This creates a new ParticleSet based on locations of all particles written
in a zarr ParticleFile at a certain time. Particle IDs are preserved if restart=True
Expand Down
2 changes: 1 addition & 1 deletion tests/test_fieldset.py
Original file line number Diff line number Diff line change
Expand Up @@ -343,7 +343,7 @@ def sampleTemp(particle, fieldset, time): # pragma: no cover
]
)

pset = ParticleSet.from_list(fieldset, pclass=MyParticle, lon=[0.5], lat=[0.5], depth=[0.5])
pset = ParticleSet(fieldset, pclass=MyParticle, lon=[0.5], lat=[0.5], depth=[0.5])
pset.execute(
AdvectionRK4_3D + pset.Kernel(sampleTemp), runtime=timedelta(hours=51), dt=timedelta(hours=dt_sign * 1)
)
Expand Down
4 changes: 2 additions & 2 deletions tests/test_particlesets.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ def test_pset_create_list_with_customvariable(fieldset):
MyParticle = Particle.add_variable("v")

v_vals = np.arange(npart)
pset = ParticleSet.from_list(fieldset, lon=lon, lat=lat, v=v_vals, pclass=MyParticle)
pset = ParticleSet(fieldset, lon=lon, lat=lat, v=v_vals, pclass=MyParticle)
assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
assert np.allclose([p.v for p in pset], v_vals, rtol=1e-12)
Expand Down Expand Up @@ -127,7 +127,7 @@ def test_pset_create_with_time(fieldset):
time = 5.0
pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle, time=time)
assert np.allclose([p.time for p in pset], time, rtol=1e-12)
pset = ParticleSet.from_list(fieldset, lon=lon, lat=lat, pclass=Particle, time=[time] * npart)
pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle, time=[time] * npart)
assert np.allclose([p.time for p in pset], time, rtol=1e-12)
pset = ParticleSet.from_line(fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=Particle, time=time)
assert np.allclose([p.time for p in pset], time, rtol=1e-12)
Expand Down
13 changes: 0 additions & 13 deletions tests/v4/test_particleset.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,19 +37,6 @@ def test_pset_create_lon_lat(fieldset):
assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)


@pytest.mark.parametrize("lonlatdepth_dtype", [np.float64, np.float32])
def test_pset_create_line(fieldset, lonlatdepth_dtype):
npart = 100
lon = np.linspace(0, 1, npart, dtype=lonlatdepth_dtype)
lat = np.linspace(1, 0, npart, dtype=lonlatdepth_dtype)
pset = ParticleSet.from_line(
fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=Particle, lonlatdepth_dtype=lonlatdepth_dtype
)
assert np.allclose([p.lon for p in pset], lon, rtol=1e-12)
assert np.allclose([p.lat for p in pset], lat, rtol=1e-12)
assert isinstance(pset[0].lat, lonlatdepth_dtype)


def test_create_empty_pset(fieldset):
pset = ParticleSet(fieldset, pclass=Particle)
assert pset.size == 0
Expand Down
2 changes: 1 addition & 1 deletion tests/v4/test_structured_gcm.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,6 @@ def test_fieldset_frompop():
# dimensions = {"lon": "lon", "lat": "lat", "time": "time"}

# fieldset = FieldSet.from_pop(filenames, variables, dimensions, mesh="flat")
# pset = ParticleSet.from_list(fieldset, Particle, lon=[3, 5, 1], lat=[3, 5, 1])
# pset = ParticleSet(fieldset, Particle, lon=[3, 5, 1], lat=[3, 5, 1])
# pset.execute(AdvectionRK4, runtime=3, dt=1)
pass
1 change: 1 addition & 0 deletions v3to4-breaking-changes.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,3 +12,4 @@ FieldSet
ParticleSet

- ParticleSet.execute() expects `numpy.datetime64`/`numpy.timedelta.64` for `runtime`, `endtime` and `dt`
- `ParticleSet.from_field()`, `ParticleSet.from_line()`, `ParticleSet.from_list()` has been removed
Loading