Preventing land 0.0 values from being sampled by a temperature sampling kernel

[TomMBIO]

Question

Question

Hi everyone,

I've been working on a dispersal simualtion that samples a temperature field in order to calculate larval growth over time. A problem I have been encountering is that, while sampling the temperature field (kernel SampleT), sometimes larvae appear to draw data from land cells (T = 0.0) resulting in erroneous dips in the temperature at the coastal boundary. The dips aren't to zero degrees however, which makes me wonder if at some point under the hood, the kernel calculates an average temperature value?

I have another kernel (Unbeaching) which forces larvae that are above land into an ocean cell, and this runs before SampleT in my kernel loop, so I assumed this sampling of land cells would not be an issue.

So in summary, I have two questions:

1) When larvae interpolate values from the field, do they take the average of the nearest grid cell points?

And,

2) If so, is there a way to prevent interpolation from grid cells with temperature values of 0.0?

I have tried a few potential solutions to this, the first was to modify the fields themselves, and assign all land cells (T= 0.0) a temperature value based on the nearest ocean cell for data to be drawn from. However this was very computationally expensive and caused a lot of memory issues. I also tried masking 0.0 values within the field to NaNs in the hopes the SampleT kernel would ignore them, but again this led to SSL errors. I think I struggle applying this logic to the many thousands of time slices my xarr dataset has.

I thought of simply adding an if statement to the SampleT kernel which ignores values of 0.0, but the issue is the sampled values are not often zero but instead a lowered number, by the zero (if that makes sense).

Any advice would be greatly appreciated. I have attached the kernel code as well as some outputs to visualise what is going on, below. Perhaps something like interp_method='linear_invdist_land_tracer' could help?

Supporting code/error messages

#---------------------- KERNELS -----------------------------

def Sample_land(particle, fieldset, time):
    """Sample if the particle is on land."""
    particle.on_land = fieldset.landmask[time, particle.depth, particle.lat, particle.lon]

def Unbeaching(particle, fieldset, time):
    """Update particle position if it is on land."""
    if particle.on_land == 1:  # If the particle is on land
        (u_land, v_land) = fieldset.UV_unbeach[time, particle.depth, particle.lat, particle.lon]  # land U and V velocities
        particle_dlon += u_land * particle.dt  # Update lon with land U velocity
        particle_dlat += v_land * particle.dt  # Update lat with land V velocity
        particle.on_land = 0  # The particle is no longer on land

def SampleT(particle, fieldset, time):
    """Sample the temperature field"""
    particle.Temp = fieldset.T[time, particle.depth, particle.lat, particle.lon]

#------------------- SIMULATE THE RELEASE OF PARTICLES -----------------------

import os

output_dir = r"_____________"

iterations = range(100) #100 larval cohorts released, each a week from the other
cohort = 0  
release_time = datetime.datetime(2023, 1, 1, 00, 00, 0)

# Track output files
output_files = []

for i in iterations:
    cohort_array = np.full(n_sites, int(cohort))  # Create cohort values
    pset = ParticleSet(
        fieldset=fieldset,
        pclass=Particle,
        lat=lats10,
        lon=lons10,
        depth=None,
        time=release_time, 
        cohort=cohort_array
    )

    output_file_name = os.path.join(output_dir, f"SA_cohort_{cohort}.zarr")
    output_files.append(output_file_name)
    output_file = pset.ParticleFile(
        name=output_file_name,
        outputdt=datetime.timedelta(minutes=60)
    )
    
    kernels = pset.Kernel(Sample_land) + pset.Kernel(Unbeaching) + \
              pset.Kernel(AdvectionRK4) + pset.Kernel(DiffusionUniformKh) + \
              pset.Kernel(SampleT) + pset.Kernel(CalculateG) + \
              pset.Kernel(GrowS) + pset.Kernel(CalculateAge) + pset.Kernel(Settlement) + \
              pset.Kernel(Kill) + pset.Kernel(CheckOutOfBounds)
    
    # Execute simulation
    pset.execute(
        kernels,
        runtime=datetime.timedelta(days=65), # Maximum achievable PLD is 64.78, so each individual will have settled before this value.
        dt=datetime.timedelta(minutes=60),
        output_file=output_file
    )
    
    # Update cohort ID and release time
    cohort += 1
    release_time += datetime.timedelta(days=7)

<img width="1140" height="717" alt="image" src="https://github.com/user-attachments/assets/55d8d609-fabf-4116-b218-a62e086d4c40" />
<img width="1140" height="717" alt="image" src="https://github.com/user-attachments/assets/c85e2825-15a9-49b6-b010-5e3267a8b04e" />

Notice the dip at the coastal boundary -- I had a look at the underlying data and this dip in 10 degrees should not exist.

Thank you for your time and consideration,

Tom

[TomMBIO]

[erikvansebille]

Hi @TomMBIO, thanks for reaching out with this question. Yes, we have support for interpolating fields that have zeros on land: it's called linear_invdist_land_tracer

See https://docs.parcels-code.org/en/latest/examples/tutorial_interpolation.html for more information

[TomMBIO]

Thanks Erik, I'll try this out and let you know how it goes.

KR
Tom

[TomMBIO]

Just an update, this worked great! (Before and After below)