data initialization

[blg0116]

Hi Angus,
Thank you for the great package! I wanted to use it for MOM6 hourly data.

I have tried to install the environment with Python3.8 but failed, then I used the Python 3.7. But on the PyPI website it's saying that the parcels package is good for python 3.8 and higher (https://pypi.org/project/parcels/#history). And I have been always getting some errors saying “NotImplementedError: Access to Field attributes are not (yet) implemented in JIT mode”. is it because the package has some inconsistency with parcels package?

Could you share the script used for the ROMS data file or do you have a version for the MOM6 data? That would be very helpful.

Sorry for my messy questions and many thanks!
Yang

[angus-g]

Hi Yang, thanks for your interest in the package. I assume you installed the library using the pip instructions? Back when I was developing this library, Parcels' only multiprocessing support was through MPI, which was completely unacceptable for what we were trying to do. It was also quite unstable around some of the features I was relying on.

For the data we were filtering, it was sufficient to use OpenMP for parallelism, for which I maintained my own fork of Parcels (at https://github.com/angus-g/parcels). Unfortunately, because I haven't looked at things for a while, this is probably quite far out of date as far as Python compatibility is concerned! I can look at updating that compatibility, it should be fairly straightforward.

As for the actual error you're hitting, I'm not quite sure what the root cause is. If you're happy to share the script you're trying to use, I might be able to provide some advice. There is an example of the kind of things to set up for variables/dimensions at https://lagrangian-filtering.readthedocs.io/en/latest/examples.html#load-roms-data. I haven't got a written-up example for MOM6, but it's largely similar.

[blg0116]

Hi Angus, thanks for your reply!
I used the Conda to install the package. Below is my script, I appreciate any suggestions you may have:

from datetime import timedelta
import os
import re

directory = '/Volumes/A4/esmg/DATA/MOM6/CAR50/exp0204/' 
pattern = r'19970101\.ocean_hourly__1997_(03[2-7])\.nc'
matching_files = []
for filename in os.listdir(directory):
    if re.match(pattern, filename):
        full_path = os.path.join(directory, filename)
        matching_files.append(full_path)

filenames = {
    "U": matching_files,
    "V": matching_files,
}
variables4 = {
    'U': 'ssu',
    'V': 'ssv',
}

dimensions2 = {
  "U":   {"lon": "xq", "lat": "yh",'time':"time"},
  "V":   {"lon": "xh", "lat": "yq",'time':"time"},
  "ssh": {"lon": "xh", "lat": "yh",'time':"time"},
}

ff = filtering.LagrangeFilter(
    'mom6simulation', filenames, variables4, dimensions2, sample_variables=["U","V"], mesh="spherical", window_size=timedelta(days=2).total_seconds(), highpass_frequency=5e-05
)
ff.make_zonally_periodic() 

with this script, I got the error

    262             funcname="sample_kernel",
    263             funcvars=["particle", "fieldset", "time"],
--> 264             funccode=f_str,
    265         )
    266 

~/anaconda3/envs/filtering/lib/python3.7/site-packages/parcels/kernel.py in __init__(self, fieldset, ptype, pyfunc, funcname, funccode, py_ast, funcvars, c_include, delete_cfiles)
    148             kernelgen = KernelGenerator(fieldset, ptype)
    149             kernel_ccode = kernelgen.generate(deepcopy(self.py_ast),
--> 150                                               self.funcvars)
    151             self.field_args = kernelgen.field_args
    152             self.vector_field_args = kernelgen.vector_field_args

~/anaconda3/envs/filtering/lib/python3.7/site-packages/parcels/codegenerator.py in generate(self, py_ast, funcvars)
    396         # Replace occurences of intrinsic objects in Python AST
    397         transformer = IntrinsicTransformer(self.fieldset, self.ptype)
--> 398         py_ast = transformer.visit(py_ast)
    399 
    400         # Untangle Pythonic tuple-assignment statements

~/anaconda3/envs/filtering/lib/python3.7/ast.py in visit(self, node)
    269         method = 'visit_' + node.__class__.__name__
    270         visitor = getattr(self, method, self.generic_visit)
--> 271         return visitor(node)
    272 
    273     def generic_visit(self, node):
NotImplementedError: Access to Field attributes are not (yet) implemented in JIT mode```

[blg0116]

Hi Angus,
I have an update. I tried one of your example (somehow I forgot where is the example is provided) as below:

%matplotlib inline
import filtering
from datetime import timedelta
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from netCDF4 import Dataset
U0=0.2
Uw=0.01
Ue=0.01
f=5e-5
k=1e-3
Le=10e3;
# Periodic grids
Lx=8*2*np.pi/k
x=np.linspace(0,Lx,200,endpoint=False)
y=np.copy(x)
ym=np.mean(y);
# Time - 2 weeks
t=np.array(range(0,(2*7*24*60**2),3600))
T,Y,X=np.meshgrid(t,y,x,indexing='ij')
# estimate of the time taken to advect 1 grid cell
(x[2]-x[1])/U0
# estimate of "frequency" of eddy
# Generate the flow, accounting for periodicity
xc=(U0*t % Lx)[:,None,None]
ut1=-2*Ue/Le*(X-xc)*np.exp(-((X-xc)**2+(Y-ym)**2)/Le**2)
vt1=2*Ue/Le*(Y-ym)*np.exp(-((X-xc)**2+(Y-ym)**2)/Le**2)
ut2=-2*Ue/Le*(X-(xc+x[-1]))*np.exp(-((X-(xc+x[-1]))**2+(Y-ym)**2)/Le**2)
vt2=2*Ue/Le*(Y-ym)*np.exp(-((X-(xc+x[-1]))**2+(Y-ym)**2)/Le**2)
ut3=-2*Ue/Le*(X+(-xc+x[-1]))*np.exp(-((X+(-xc+x[-1]))**2+(Y-ym)**2)/Le**2)
vt3=2*Ue/Le*(Y-ym)*np.exp(-((X+(-xc+x[-1]))**2+(Y-ym)**2)/Le**2)
ue=(ut1+ut2+ut3)
ve=(vt1+vt2+vt3)
# total flow
U=U0+Uw*np.sin(k*X)+ue
V=f*Uw/(k*U0)*np.cos(k*X)+ve
ti=100;
plt.subplot(2,1,1)
plt.pcolor(x,y,U[ti,:,:])
plt.title('u')
plt.colorbar()
plt.subplot(2,1,2)
plt.pcolor(x,y,V[ti,:,:])
plt.title('v');
plt.colorbar()
ncfile_name="analytic_example_input.nc"
ncfile = Dataset(ncfile_name, 'w', format='NETCDF4')
ncfile.description = 'Synthetic data for Lagrangian Filtering test'
# dimensions
ncfile.createDimension('time', None)
ncfile.createDimension('x', np.size(x))
ncfile.createDimension('y', np.size(y))
# variables
time = ncfile.createVariable('time', 'f8', ('time',))
x1 = ncfile.createVariable('x', 'f8', ('x',))
y1 = ncfile.createVariable('y', 'f8', ('y',))
U1 = ncfile.createVariable('u', 'f8', ('time', 'y', 'x'))
V1 = ncfile.createVariable('v', 'f8', ('time', 'y', 'x'))
U1[:] = U
V1[:] = V
time[:] = t
x1[:] = x
y1[:] = y
ncfile.close()

filenames = {
"U": ncfile_name,
"V": ncfile_name
}
variables = {"U": "u", "V": "v"}
dimensions = {"lon": "x", "lat": "y", "time": "time"}
ff = filtering.LagrangeFilter(
"analytic_example_out1", filenames, variables, dimensions,
sample_variables=["U","V"], mesh="flat",
window_size=timedelta(days=3.5).total_seconds(),
highpass_frequency=f
)
# periodic domain
ff.make_zonally_periodic()
ff.make_meridionally_periodic()
ff.advection_dt=600
t1=timedelta(days=7).total_seconds()
ff.filter(times=[t1])

I got the same error 'NotImplementedError: Access to Field attributes are not (yet) implemented in JIT mode'
I used the ' conda create -n filtering -c angus-g -c conda-forge lagrangian-filtering' to install the package and the installed parcels version is 'parcels 2.2.1.dev63+g9e8533e3 '. I guess it's the compatibility issue as you suggested.