DOCS: Tutorial on how to get repeated releases now that repeatdt is removed

[VeckoTheGecko]

This can involve both the nested for loop approach as well as some numpy magic to get the same goal :)

[michaeldenes]

As a comment, I think highlighting even more to users that particles releases are actually a set of (longitude, latitude, depth, time) coordinates would be beneficial (and not just simple (longitude,latitude) pairs where parcels assumes depth=0 and time=0). (and then of course showing tips on how to construct lists of each of these).

[nilodna]

I’m sharing the code I’m using to seed particles on a vertical section (25 particles arranged in a 5x5 matrix). Because of the section’s orientation, things get a bit trickier, so the code still needs some revision to make it more general. Hopefully it might be useful.

import numpy as np

ndays = 1 # should be the number of timesteps you want to seed

lon_bounds = [-49.5, -45]
lat_bounds = [0,4]
dep_bounds = [5,100]

lat = np.linspace(*lat_bounds, 5)
lon = np.linspace(*lon_bounds, 5)
dep = np.linspace(*dep_bounds, 5)
dates = np.arange(ndays) * 86400  # seconds, to seed every day

# build grids (ndays, dep, lon)
lon_grid, dep_grid, time_grid = np.meshgrid(lon, dep, dates, indexing="ij")

# match shapes
lons   = lon_grid.transpose(2,1,0)
deps   = dep_grid.transpose(2,1,0)
times  = time_grid.transpose(2,1,0)

# repeating the list of latitudes to match the number of depths
lats   = np.broadcast_to(lat[np.newaxis, np.newaxis, :], lons.shape)

# then flat all matrices for Parcels and save as a npy file
modellons = lons.flatten()
modellats = lats.flatten()
modeldepths = deps.flatten()
modeltimes = times.flatten()

np.savez("seeding_coords.npz", 
         lons=modellons, 
         lats=modellats, 
         deps=modeldepths, 
         time=modeltimes)

Then to illustrate:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # needed for 3D projection

# Remove the extra axis (shape is (1,5,5), we want (5,5))
lons, lats, deps = lons[0], lats[0], deps[0]

# Flatten to 1D arrays for scatter
x = lons.flatten()
y = lats.flatten()
z = deps.flatten()

# nt = number of timesteps just to show how the matrices are the same
nt = 4
fig, axes = plt.subplots(1, nt, subplot_kw={'projection':'3d'}, figsize=(14,5))
for i,ax in enumerate(axes.ravel()):
    
    # Scatter plot
    ax.scatter(x, y, z, c=z, cmap="viridis", s=50, edgecolor="k")
    ax.set_title(f"t = {i*86400}")
    if i == 3:
        ax.set_zlabel("Depth")
        
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    
    ax.invert_zaxis()

plt.show()

[reint-fischer]

Implemented in #2357

We have also now updated many tutorials to include initial values for time and z as opposed to just lat and lon.