Trouble animating my parcels output

[TomMBIO]

Hi there,

I am trying to animate my parcels output but am not having any luck. I would like to animate each time step within the zarray but am not sure how. I think part of the problem is that I'm calling to multidimensional data when the plotting packages require one dimensional data but I'm not sure at all.

I've attached my code, the model runs fine but all the plotting things can be found at the bottom. Can you shed some light on this? I'm a bit lost.

(By the way, what is the best format to upload my code in for future reference?)

Thanks,
Tom

Running the model

import math
from datetime import timedelta
from operator import attrgetter

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import trajan as ta
import xarray as xr
from IPython.display import HTML
from matplotlib.animation import FuncAnimation

from parcels import (
AdvectionRK4,
DiffusionUniformKh,
FieldSet,
JITParticle,
ParticleSet,
Variable,
)

import netCDF4 as nc
import numpy as np
import itertools as it

Open the NetCDF file in read mode

filename = 'C:/Users/benth/Documents/PhD/NetCDF/vector_data.nc'

try:
nc_data = nc.Dataset(filename, 'r')
print("File opened successfully.")
except Exception as e:
print("Error opening file:", e)

#this section will show you whats in the netcdf file
#-------------------------------------------------

Print global attributes

print("Global attributes:")
for attr in nc_data.ncattrs():
print(f"{attr}: {getattr(nc_data, attr)}")

Print dimensions

print("\nDimensions:")
print(list(nc_data.dimensions.keys()))

Print variables

print("\nVariables:")
print(list(nc_data.variables.keys()))

Close the NetCDF file

nc_data.close()

#----------------------------------------------

filenames = {
"U": filename ,
"V": filename ,
}

variables = {"U": "uo", "V": "vo"}

dimensions = {"lat": "latitude", "lon": "longitude", "time": "time"}

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)

Add diffusivity. Testing tiny diffusivity to prevent from advecting out of bounds.

kh = 10
fieldset.add_constant_field("Kh_zonal", kh, mesh = "spherical")
fieldset.add_constant_field('Kh_meridional', kh, mesh='spherical')

#creating list of lat and longs for pset release points. These were extracted from rockyshore grid made in QGIS.
df = pd.read_csv(r"C:\Users\benth\Documents\PhD\GIS\UK littorals\rocky_gridpoints_pruned.csv")
particle_lat = df['Latitude'].tolist() #to extract the lats and longs specifically
particle_lon = df['Longitude'].tolist()

Define a custom kernel to update the particle's time attribute

def update_time(particle, fieldset, time):
particle.time = time

pset = ParticleSet.from_list(
fieldset=fieldset,
pclass=MyParticle,
lat=particle_lat,
lon=particle_lon,
depth=None,
)

output_file = pset.ParticleFile( name=r"C:\Users\benth\Documents\PhD\ocean_parcels\parcels_testing\Particles60dt960day.zarr", outputdt=timedelta(minutes=5))

Define kernels, which are looped over and applied to each advecting particle

kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(DiffusionUniformKh)

#Run particles
pset.execute(
kernels,
runtime=timedelta(hours=960),
dt=timedelta(minutes=60),
output_file=output_file)

ds = xr.open_zarr(r"C:\Users\benth\Documents\PhD\ocean_parcels\parcels_testing\Particles60dt960day.zarr")
ds.traj.plot(margin=0.4)
plt.show()

from IPython.display import HTML
from matplotlib.animation import FuncAnimation

Now to animate the output

%%capture
fig = plt.figure(figsize=(5, 5), constrained_layout=True)
ax = fig.add_subplot()
ds = xr.open_zarr(r"C:\Users\benth\Documents\PhD\ocean_parcels\parcels_testing\Particles60dt960day.zarr") #particleset output file

outputdt = timedelta(hours=1)

timerange in nanoseconds

timerange = np.arange(
np.nanmin(ds.time.values),
np.nanmax(ds.time.values) + np.timedelta64(outputdt),
outputdt,
)

ax.set_ylabel("Meridional distance [m]")
ax.set_xlabel("Zonal distance [m]")
ax.set_xlim(min(ds.lat.values.flatten()), max(ds.lat.values.flatten()))
ax.set_ylim(min(ds.lon.values.flatten()), max(ds.lon.values.flatten())) #no - axes
plt.xticks(rotation=45) #changes orientation of x axis data to be 45*

Indices of the data where time = 0

time_id = np.where(ds.time.values == timerange[0])

scatter = ax.scatter(
ds.lon.values[time_id], ds.lat.values[time_id]
)

t = str(timerange[0].astype("timedelta64[h]"))
title = ax.set_title("Particles at t = " + t)

def animate(i):
t = str(timerange[i].astype("timedelta64[h]"))
title.set_text("Particles at t = " + t)

time_id = np.where(ds.time.values == timerange[i])
scatter.set_offsets(
    np.c_[ds.lon.values[time_id], ds.lat.values[time_id]]
)

anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)

HTML(anim.to_jshtml())

[michaeldenes]

Hi @TomMBIO, what exactly are you trying to animate? Just particles updating to their new positions at each timestep? We have a few examples in the OceanParcels gallery here: https://oceanparcels.org/#gallery as well as in the homepage animation tutorial here: https://docs.oceanparcels.org/en/latest/examples/documentation_homepage_animation.html and the particle-particle interaction tutorial here: https://docs.oceanparcels.org/en/latest/examples/tutorial_interaction.html

Let me know how you go, cheers,
Michael

[TomMBIO]

This looks promising, I'll let you know how I get on!

And yes, just trying to animate the particles moving through each timestep.
Cheers,
Tom