Correct sampling of constant 2D-field in time-varying 3D-run

[RomanIMMSP]

Dear Parcels community,

I met a problem, which probably has a simple solution.
I want to track the 2D-field (ocean depth in my case) along the particle trajectory.
I define the H (depth) variable as

filenames = { ...
             'H': {'lon': bath_file, 'lat': bath_file, 'data': bath_file}  }
variables = {..., 'H': 'h'}
dimensions = {...
              'H': {'lon': 'lon_rho', 'lat': 'lat_rho'}  }

So H does not change in time.

How correctly can I sample the H along the particle trajectory? I try to do this in a such a way

def SampleH(particle, fieldset, time):
     particle.h = fieldset.H[0, 0, particle.lat, particle.lon]

But the model gives an error time origins of different grids must be have the same time
As I understand, this is because dimensions of H do not include time.
I don't know if it is important, but this error appears when time_periodic option is ON.

In this case one more issue.
The example in the Field sampling tutorial gives
particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]
The example in the Tutorial on Parcels tracer interpolation methods gives that
particle.p = fieldset.P[time+particle.dt, particle.depth, particle.lat, particle.lon]
As I understand, example with time+particle.dt is more correct.

Thanks in advance to everyone who will react to this message.

[erikvansebille]

Hi @RomanIMMSP, thanks for reporting. It could indeed very well be that you run into this problem because you use time_periodic=True. That option is incompatible with the allow-Time_extrapolation=True that is needed for sampling time-constant fields. Could you confirm that the error disappears if you set periodic time to False?

[RomanIMMSP]

Hi @erikvansebille. Yes, the error disappears when time_periodic is False.

[erikvansebille]

OK, this is a use-case we haven't seen before. I will see if I can create a hotfix to support sampling time-constant fields when using time periodic FieldSet. But it may be a few days until I find the time for it...

[erikvansebille]

I looked again at your issue, but I can't reproduce it on a simple unit test script, see 2555163. So perhaps the problem has to do more with the time_origins than with the time_extrapolation?

[RomanIMMSP]

I have solved the problem by adding the dimension ocean_time in the H variable.

  time=ds.ocean_time.shape[0]
  hh=ds.h.expand_dims({'ocean_time': time}).assign_coords({'ocean_time': ds.ocean_time})

After that, I have the same array H at each time step. And the model works without errors.

Thanks for the support.

[LAWTIDEBOY]

Dear @erikvansebille, all,

I am sorry but I am also struggling quite a bit with the sampling of a constant field such as a bathymetry field. Maybe I have done stupid mistakes (see enclosed my notebook and data for better reproducibility: https://we.tl/t-PnmssqFQO5).

The script "works" but the sampling field somehow is not consistent with the sampling field.
In other words, the depths along the trajectories are not matching the actual depths after visual inspection.

I already ran into the documentation and this forum, tried almost everything (getting rid of other kernels, disallow time extrapolation, etc) but I arrive at a point where I would need a little help. That would be lovely (I would like to avoid a time-varying transformation of the bathymetric data as much as possible as proposed by Roman).

Let me describe my workflow.
I basically create the Field Set as followed:

files = sorted(glob(f"ARC_uo_glor_*_01.nc"))
bathy_file = 'etopo025.nc'

filenames = {'U': files,
             'V': files,
             'H': bathy_file}

variables = {'U': 'uo_glor',
             'V': 'vo_glor',
             'H': 'Z'}

full_dims = {'lat': 'latitude',
              'lon': 'longitude',
              'time': 'time'}

dimensions = {'U': full_dims, 'V': full_dims, 'H': {'lat': 'latitude', 'lon': 'longitude'}}

# Create patch of initial particles

[xi,yi] = np.meshgrid(np.arange(-20,11,1),np.arange(79,80.25,0.25), indexing='ij')
y = np.random.randn(xi.shape[0],xi.shape[1])

xi = xi.reshape(-1, 1)
yi = yi.reshape(-1, 1)

# Create Field Set
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True)
fieldset.add_constant('halo_west', fieldset.U.grid.lon[0])
fieldset.add_constant('halo_east', fieldset.U.grid.lon[-1])
fieldset.add_periodic_halo(zonal=True)

Then I create the sampling Kernel:

# Create Kernel for Sampling constant bathymetry field

class SampleParticle(JITParticle):            # Define a new particle class
    H = Variable('h', dtype=np.float32, initial=0.)  # Variable 'depth' 

pset = ParticleSet(
    fieldset=fieldset, 
    pclass=SampleParticle, 
    lon=xi, # a vector of release longitudes      
    lat=yi, # a vector of release latitudes
    time=datetime(2000, 8, 1)
)

def SampleH(particle, fieldset, time):
    particle.h = fieldset.H[0,0,particle.lat, particle.lon]

And finally, I derive:

# Derive Lagrangian trajectories

output_file = pset.ParticleFile(name="Particles_2000.zarr",outputdt=timedelta(days=1))

pset.execute([AdvectionRK4, SampleH, periodicBC, DeleteErrorParticle],
             runtime=timedelta(days=28),
             dt=-timedelta(days=1),# the timestep of the kernel
             output_file=output_file) #

As you can see below, the bathymetry and the depth sampled along the trajectory are not matching:

[erikvansebille]

Thanks for reporting, @RomanIMMSP. How strange that this doesn't work!

A few questions/comments:

1. I see you use H = Variable('h', ... in your custom particle class. Better to use the same symbol for the variable name and the assignment, so both lowercase? Although I don't think this makes the difference
2. What if you remove the periodicBC kernel from the pset.execute()? Same problem then?
3. What if you don't do a pset.execute() at all, but simply type print(fieldset.H[0, 0, 79, -5]). Does that give you the approx -300 value that suggests in the plot? And what if you type print(fieldset.H[0, 0, 79, -18])?

[LAWTIDEBOY]

Dear @erikvansebille,

Thank you so much for your swift answer.

The two first options I already tried before without effect.

The third one however is indeed showing wrong outputs. I found:

print(fieldset.H[0, 0, 79, -5]) = -49
print(fieldset.H[0, 0, 79, -18]) = 0

The first one should print a value associated with the basin, i.e. thousands of meters deep.
I am out of clue unfortunatly.

PS: I must say that for the sake of my project, I finally used @RomanIMMSP s solution, with a "fake" time-varying array of bathymetry (i.e. a repeated bathymetry) and this works nicely as a classic sampling although it is wasting storage for long advections.

[erikvansebille]

Hmm, and what if you treat it as a 3D array? So

print(fieldset.H[0, 79, -5])
print(fieldset.H[0, 79, -18])

[LAWTIDEBOY]

In this case I have an error saying there is one missing argument :/

TypeError: Field.eval() missing 1 required positional argument: 'x'

[erikvansebille]

Ah ok, that’s clear. And what if you use the fix suggested by @RomanIMMSP above, to set allow_time_extrapolation=False? That seemed to work in his case?

[LAWTIDEBOY]

Hello, thank again for the support !

I have tried allow_time_extrapolation=False alone and did not work.

For a short advection such as in this example (30days) this is quite straightforward indeed to work out with Roman 's solution combined with allow_time_extrapolation=False. It works of course.

However, I do 5 to 6 years advection experiments, with daily time steps.

I found it easier to keep working "from netcdf" in this case. So I just created new source netcdf files that include the velocities together with a repeated (in time) bathymetry, and the same date time index. In this way, although it wastes some storage, it works very nicely. This is the only successful work around I have found in my case.

So if I understand well, there is two solutions for now:

• either you work from datastructure like Roman and you add a time dimension to the constant field
• either you work fromn netcdf files like myself, and artificially modify the source netcdf to add the repeated (in time) constant field.