T and S output with parcels in POP

[norafried]

Hi,
I'm trying to get T and S properties along the particles tracks with POP data. I followed this tutorial and got until the kernel. Then got an error concerning the grid. I'm assuming the Bgrid in POP might be the reason for that. Has tried to do that before in POP and can direct me to an example?

Thanks a lot!

[erikvansebille]

Dear @norafried, thanks for askign this question. But for me to best help you, I'd need to know a bit more about your case. Which code exactly did you run? What was the full error? I can't help you if you don't provide us with this more detailed information

[norafried]

Hi @erikvansebille, thanks for your reply and yes I'm sorry that I didn't make that any clearer. I followed this tutorial:

https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_sampling.ipynb?flush_cache=true

I get the error in box 5 of the tutorial when running the particle set with the kernel in the very last line.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TypeError Traceback (most recent call last)
Cell In[13], line 3
1 # define all kernels to be executed on particles
2 # kernels = pset.Kernel(AdvectionRK4_3D) #+ pset.Kernel(periodicBC)
----> 3 kernels = pset.Kernel(SampleUVTS)

File ~/miniconda3/envs/parcels_215/lib/python3.11/site-packages/parcels/particleset.py:763, in ParticleSet.Kernel(self, pyfunc, c_include, delete_cfiles)
758 def Kernel(self, pyfunc, c_include="", delete_cfiles=True):
759 """Wrapper method to convert a pyfunc into a :class:parcels.kernel.Kernel object
760 based on fieldset and ptype of the ParticleSet
761 :param delete_cfiles: Boolean whether to delete the C-files after compilation in JIT mode (default is True)
762 """
--> 763 return Kernel(self.fieldset, self.ptype, pyfunc=pyfunc, c_include=c_include,
764 delete_cfiles=delete_cfiles)

File ~/miniconda3/envs/parcels_215/lib/python3.11/site-packages/parcels/kernel.py:134, in Kernel.init(self, fieldset, ptype, pyfunc, funcname, funccode, py_ast, funcvars, c_include, delete_cfiles)
132 if self.ptype.uses_jit:
133 kernelgen = KernelGenerator(fieldset, ptype)
--> 134 kernel_ccode = kernelgen.generate(deepcopy(self.py_ast),
135 self.funcvars)
136 self.field_args = kernelgen.field_args
137 self.vector_field_args = kernelgen.vector_field_args

File ~/miniconda3/envs/parcels_215/lib/python3.11/site-packages/parcels/codegenerator.py:391, in KernelGenerator.generate(self, py_ast, funcvars)
388 self.docstr = ast.get_docstring(py_ast)
390 # Generate C-code for all nodes in the Python AST
--> 391 self.visit(py_ast)
392 self.ccode = py_ast.ccode
394 # Insert variable declarations for non-instrinsics
395 # Make sure that repeated variables are not declared more than
396 # once. If variables occur in multiple Kernels, give a warning

File ~/miniconda3/envs/parcels_215/lib/python3.11/ast.py:410, in NodeVisitor.visit(self, node)
408 method = 'visit_' + node.class.name
409 visitor = getattr(self, method, self.generic_visit)
--> 410 return visitor(node)

File ~/miniconda3/envs/parcels_215/lib/python3.11/site-packages/parcels/codegenerator.py:422, in KernelGenerator.visit_FunctionDef(self, node)
420 for stmt in node.body:
421 if not (hasattr(stmt, 'value') and type(stmt.value) is ast.Str): # ignore docstrings
--> 422 self.visit(stmt)
424 # Create function declaration and argument list
425 decl = c.Static(c.DeclSpecifier(c.Value("ErrorCode", node.name), spec='inline'))

File ~/miniconda3/envs/parcels_215/lib/python3.11/ast.py:410, in NodeVisitor.visit(self, node)
408 method = 'visit_' + node.class.name
409 visitor = getattr(self, method, self.generic_visit)
--> 410 return visitor(node)

File ~/miniconda3/envs/parcels_215/lib/python3.11/site-packages/parcels/codegenerator.py:747, in KernelGenerator.visit_VectorFieldEvalNode(self, node)
745 self.visit(node.field)
746 self.visit(node.args)
--> 747 ccode_eval = node.field.obj.ccode_eval(node.var, node.var2, node.var3,
748 node.field.obj.U, node.field.obj.V, node.field.obj.W,
749 *node.args.ccode)
750 if node.field.obj.U.interp_method != 'cgrid_velocity':
751 ccode_conv1 = node.field.obj.U.ccode_convert(*node.args.ccode)

TypeError: VectorField.ccode_eval() takes 11 positional arguments but 16 were given

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Let me know in case you need more info.

[erikvansebille]

Hmm, the error says that the problem is in your SampleUVTS Kernel. Can you provide the code for that?

[erikvansebille]

OK, then I understand the error. This Sampling Kernel is incompatible with version 2.1.5. You would need to do something like

def SampleUVTS(particle,fieldset,time):
    particle.u = fieldset.U[time,particle.depth,particle.lat,particle.lon]
    particle.v = fieldset.V[time,particle.depth,particle.lat,particle.lon]
    particle.t = fieldset.T[time,particle.depth,particle.lat,particle.lon]
    particle.s = fieldset.S[time,particle.depth,particle.lat,particle.lon]

[norafried]

This made the Kernel 'work'. But if I proceed further and inititalize the field at dt=0 the kernel dies. It doesn't give me an error, it just dies.

pset.execute(pset.Kernel(SampleUVTS),dt=0)

If I don't initialize the field then I do get an output for T and S but just containing zeros. I guess one small thing is still incorrect?! Do you have an idea what I could check?

[erikvansebille]

It could very well be that dt=0 was not supported yet in Parcels v2.1.5? If you want to check that out, you could go through the release notes of all the versions since then.

In any case, you are running a two-year old version of the code. And I'm sorry, but we don't have the people-power to keep supporting old versions, so I can't really help you except to suggest that you upgrade to the latest version.

[norafried]

I tried running the script with Parcels 2.4.0. The error I'm getting now is similar to the one I got in the very beginning and which made me switch to an earlier Parcels version. When I'm executing the particle set it always returns me:

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 43 and the array at index 1 has size 42

Which I remember being a problem related to POP. Is there a way how to solve this in the new Parcels version? Why did this work in the earlier releases?

Thanks for your help.

[erikvansebille]

I wouldn't know how to solve your problem without seeing your code. Perhaps send me the code and a bit of the data, and I will have a look? We have some people here working with POP data, so Parcels should work for that output too