-
QuestionQuestionI'm trying to set up Pacels with a combination of field sets from different sources (and different grids). The sea surface currents come from CMEMS (https://data.marine.copernicus.eu/product/NWSHELF_MULTIYEAR_PHY_004_009/description) and wind velocities come from ERA5 (https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels?tab=overview). With the code below, I got a segmentation fault when adding the two ERA5 wind velocity fields. Surprisingly, I don't have the problem if I only add one of the two. Even more surprisingly, I don't have the segmentation fault every time (two successful runs out of dozens of attempts). (I have successful setups combining the sea surface currents dataset from CMEMS with Stokes Drift velocities from another CMEMS product. I did not include that here, to keep the example as small as possible.) Input files to reproduce the error:
I would really appreciate if someone could help me figure out what my problem is. Supporting code/error messagesfrom datetime import datetime, timedelta
from glob import glob
import numpy as np
import parcels
###############
# Field sets. #
###############
# Import sea surface currents (CMEMS).
filenames = {'U': 'ns_ssc_2007_09.nc', 'V': 'ns_ssc_2007_09.nc'}
variables = {'U': 'uo', 'V': 'vo'}
dimensions = {'time': 'time', 'lat': 'latitude', 'lon': 'longitude'}
fieldset = parcels.FieldSet.from_netcdf(filenames, variables, dimensions)
# Import wind velocities (ERA5).
filenames = 'ns_wind_2007_09.nc'
variables = {'U_W': 'u10', 'V_W': 'v10'}
dimensions = {'time': 'valid_time', 'lat': 'latitude', 'lon': 'longitude'}
fieldset_W = parcels.FieldSet.from_netcdf(filenames, variables, dimensions)
fieldset.add_field(fieldset_W.U_W, name = 'U_W')
fieldset.add_field(fieldset_W.V_W, name = 'V_W')
#################
# Particle Set. #
#################
# Define release position.
lat = [49.85] * 24
lon = [-2.82] * 24
# Define release date.
date_init_fieldset = datetime(2007, 9, 1)
date_release = datetime(2007, 9, 13)
# Release every hour of the release date.
time = (date_release - date_init_fieldset).total_seconds() + np.arange(24) * timedelta(hours = 1).total_seconds()
# Create particle set.
pset = parcels.ParticleSet(
fieldset = fieldset,
pclass = parcels.JITParticle,
lon = lon,
lat = lat,
time = time
)
############
# Kernels. #
############
# Wind drift.
def WindDrift(particle, fieldset, time):
# Windage coefficient.
alpha = 0.01
# To convert from m/s to deg/s.
R_earth = 6371000.0
deg_to_m_lat = math.pi * R_earth / 180.0
lat_rad = particle.lat * math.pi / 180.0
deg_to_m_lon = deg_to_m_lat * math.cos(lat_rad)
# Particle positions are updated only after evaluating all terms.
particle_dlon += alpha * fieldset.U_W[particle] / deg_to_m_lon * particle.dt
particle_dlat += alpha * fieldset.V_W[particle] / deg_to_m_lat * particle.dt
# Create kernels.
kernels = pset.Kernel([parcels.AdvectionRK4, WindDrift])
########
# Run. #
########
# Output file.
output_file = pset.ParticleFile('out.zarr', outputdt = timedelta(hours = 1))
# Execution.ParcelsRandom
date_final_fieldset = datetime(2007, 9, 29)
runtime = (date_final_fieldset - date_release).total_seconds()
pset.execute(
kernels,
runtime = runtime,
dt = timedelta(minutes = 10),
output_file = output_file
)This is the error I got (I'm running Parcels v3.1.4 on macOS): |
Beta Was this translation helpful? Give feedback.
Answered by
ogourgue
Oct 3, 2025
Replies: 4 comments
-
|
For context, I was discussing this with @ogourgue in person during the conference and I didn't know what the problem was. Just confirming @ogourgue , which version of parcels are you using? ( |
Beta Was this translation helpful? Give feedback.
0 replies
-
|
This is the version of Parcels: This is the entire environment export: |
Beta Was this translation helpful? Give feedback.
0 replies
-
|
It seems that I solved my issue by first addressing the warning appearing before the segmentation fault: I flipped the latitudes of the entire ERA5 dataset in preprocessing, and both the warning and the segmentation fault disappeared. I keep the discussion open in case you think that there might remain a bug in Parcels, but my problem is solved. Thanks again to the entire team for the very interesting workshop! I come back home with many ideas and a lot of enthusiasm for Lagrangian analyses and modeling :-) |
Beta Was this translation helpful? Give feedback.
0 replies
Answer selected by
ogourgue
-
|
Thanks @ogourgue for reporting. Glad you found the root of the issue. This will not be an issue in v4 anymore, as the searching of the grids is much more robust in the new version. But while you're on v3; be careful when combining grids with different directions for latitude! |
Beta Was this translation helpful? Give feedback.
0 replies
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
It seems that I solved my issue by first addressing the warning appearing before the segmentation fault:
I flipped the latitudes of the entire ERA5 dataset in preprocessing, and both the warning and the segmentation fault disappeared.
I keep the discussion open in case you think that there might remain a bug in Parcels, but my problem is solved.
Thanks again to the entire team for the very interesting workshop! I come back home with many ideas and a lot of enthusiasm for Lagrangian analyses and modeling :-)