Noncontinuous grids

[hblumberg]

My team has noticed that if we construct a Fieldset using NetCDF data in which longitudes are in the range [0, 360), we see out of bounds errors when a particle tries to move across the Prime Meridian. This is the case even when the longitudes are ordered "correctly" (e.g. [358, 359, 0, 1, 2]).

We could change our longitudes to be in the range [-180, 180), but then we'd have the problem along the International Date Line.

Our questions are:

1. Is this the intended behavior for noncontinous grids?
2. If so, how do folks usually work around this when working with longitudes?

Thanks in advance!

[erikvansebille]

Thanks for your question, @hblumberg. Indeed, Parcels expects longitudes to be monotonously increasing. In your example about, you would use [358, 359, 360, 361, 362]. This can be done by changing the longitude vector before (or after) the FieldSet creation, but certainly before the particles execution. Does that work?

Note that for dealing with periodic boundaries, you may then also need to add a halo. See for example the tutorial at https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_periodic_boundaries.ipynb

An alternative to using a halo is to change the kernels themselves, and add modulo statements to for example a custom Advection Kernel.

[hblumberg]

Thank you so much for your quick response (and I’m sorry for the delay on my end)!

Yes, modifying the FieldSet to have increasing longitudes works!

Creating a custom kernel would be an alternative to modifying the FieldSet -- is that correct? If so, is this the right idea for a custom version of AdvectionRK4?

def AdvectionRK4Modified(particle, fieldset, time):
    (u1, v1) = fieldset.UV[particle]
    lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt)

    if lon1 < 0:
        lon1 += 360
    if lon1 >= 360:
        lon1 -= 360

    (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle]
    lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt)

    if lon2 < 0:
        lon2 += 360
    if lon2 >= 360:
        lon2 -= 360

    (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle]
    lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt)

    if lon3 < 0:
        lon3 += 360
    if lon3 >= 360:
        lon3 -= 360

    (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
    particle.lon += (u1 + 2*u2 + 2*u3 + u4) / 6. * particle.dt
    particle.lat += (v1 + 2*v2 + 2*v3 + v4) / 6. * particle.dt

    if particle.lon < 0:
        particle.lon += 360
    if particle.lon >= 360:
        particle.lon -= 360

This doesn’t seem to work for me, but I'm thinking it might be because my longitudes are really […, 359.84, 359.92, 0, 0.08, 0.16, …]. So maybe the particles are getting lost between 359.92 and 0? Just wanted to check my thinking on the approach before I keep going in this direction.

Thanks again for your help!!

[erikvansebille]

I don't think this kernel will seamlessly work, because the problem is more under the hood, in finding which grid cell a particle is in. If the longitudes are not monotonically increasing, this builtin-function may break. So much safer to fix your fields longitudes, than to change your Kernels.

[hblumberg]

Got it, thank you!

[vbalza]

Hi Erik! Hannah and I implemented your suggestion to ensure field longitudes are monotonically increasing. After kernel execution, we want to convert the longitude values of the output particle trajectories back to ensure their range is [0, 360).

However, we are not sure how to access and modify output longitudes before exporting the file. Currently, we reload the file as an xarray dataset, modify the longitudes, and then overwrite the original output file:

    if lon_max < lon_min:
        with xr.open_dataset(deployment.output_filename) as ds:
            data = ds.load()
        data['lon'].values = np.where((data['lon'].values >= 360.0),
                                      data['lon'].values - 360.0,
                                      data['lon'].values)
        data.to_netcdf(deployment.output_filename, mode='w')
        data.close()

While this execution runs, we wanted to see if you had any recommendations regarding how to go about modifying output values to avoid reloading/overwriting the original output file. Thanks in advance for your help!

[erikvansebille]

Hi @vbalza and @hblumberg. Not sure how much advice I can give on working with output files. Of course, you can keep them in the original Parcels format and change the longitudes using your np.where() code every time you read the files in for analysis. Or you could overwrite, like you suggest above. There is no (simple) way to change values in Parcels files themselves

[vbalza]

Got it! Thank you very much for the quick response.

[vbalza]

Hi Erik,

@hblumberg and I have a follow-up question regarding fieldset modification. First, we verified that when we adjust the dataset directly (before fieldset creation) by loading the data into an xarray and adjusting the longitude coordinates to be monotonically increasing, the simulation works as intended. However, if we instead load the original dataset in which the longitude bounds are, for example, [350, 10], and then try to adjust the fieldset after fieldset creation, the simulation breaks.

Please find this approach below for a dataset with bounds [350, 10]:

# fieldset has two fields, U and V.
lon_min = fieldset.U.lon[0] # left longitude bound, 350
lon_max = fieldset.U.lon[-1] # right longitude bound, 10

if lon_max < lon_min:
    # Prime Meridian crossed
    # Convert Field objects
    fieldset.U.lon = np.where(fieldset.U.lon < lon_min, 
                              fieldset.U.lon + 360.0, 
                              fieldset.U.lon)
    fieldset.V.lon = np.where(fieldset.V.lon < lon_min, 
                              fieldset.V.lon + 360.0, 
                              fieldset.V.lon)

    # Convert Grid
    fieldset.U.grid.lon = np.where(fieldset.U.grid.lon < lon_min, 
                                   fieldset.U.grid.lon + 360.0, 
                                   fieldset.U.grid.lon)
    fieldset.V.grid.lon = np.where(fieldset.V.grid.lon < lon_min, 
                                   fieldset.V.grid.lon + 360.0, 
                                   fieldset.V.grid.lon)

# some other code [...]

pset = ParticleSet.from_list(fieldset=fieldset, 
                             lon=deployment.loc_policy.lon_E_list,  # release longitudes (where the longitudes are monotonically increasing)
                             # <other parameters>)

# some other code [...]

pset.execute(
    # <code executing kernel>)

Does parcels store the longitude coordinates elsewhere in the fieldset? Are there other attributes we need to be modifying beyond fieldset.U.lon, fieldset.V.lon, fieldset.U.grid.lon, and fieldset.V.grid.lon?

We are happy to share a minimal, executable example if that would be helpful in answering our questions. Thank you in advance for any clarification!

[erikvansebille]

What is the error that you're getting when you use this code? And if you can indeed share a minimal example, then I'll be able to try if I can fix it

[vbalza]

Hi Erik,

Thanks for getting back to us. There is no error message, as the simulation when we manually modify the fieldset runs. However, when we check the output for a dataset with longitude bounds [350, 10], for example, it seems particles to the right of the Prime Meridian (with a starting point between 0 and 10) do not simulate—even after changing the fieldset longitudes to be monotonically increasing.

I am attaching a simple example notebook parcels_test.ipynb that uses HYCOM water velocity data test.nc visualizing the above:
parcels_test.zip

Please let us know if you have any questions.

Thank you very much!
Valeria

[erikvansebille]

Thanks for this minimal example; it helped me to very efficiently find the issue/bug!

Fortunately, it's quite simple to solve in your script. After you changed the Lons and lats in the FieldSet, you only need to add

fieldset2.U.grid.lonlat_minmax[0] = min(fieldset2.U.grid.lon)
fieldset2.U.grid.lonlat_minmax[1] = max(fieldset2.U.grid.lon)

Then, it should work!

[vbalza]

Thank you so much for the quick response, Erik! That works.