dealing with particles on land

[mollyjames2]

I am trying to write a function to acknowledge if a particle is on land. My end goal is to create a csv of the date/time that a particle hits land and its lat/lon.

so far i have created a land mask that matches the dimensions of the U,V fields as follows:

# Open the files containing the U and V velocities
ufields = fieldset.U.data[0,0,:,:]
vfields = fieldset.V.data[0,0,:,:]

landmask = np.logical_or(np.ma.masked_equal(ufields,0.).mask, # Creates a masked array
                         np.ma.masked_equal(vfields,0.).mask)


landmask = landmask.astype('int') # Convert to integers


landmask = landmask[ np.newaxis, :, :]

# Repeat array1 along the new axes to match the desired shape (31, 424, 324)
landmask = np.tile(landmask, (31, 1, 1))



fieldset.add_field(Field('landmask', landmask, lon=fieldset.U.grid.lon, lat=fieldset.U.grid.lat, depth=fieldset.U.grid.depth, mesh='spherical')) # Add to fieldset.

This landmask looks fine (see image)

I have defined a function as follows:

def Sample_land(particle, fieldset, time):
    particle.on_land = fieldset.landmask[time, particle.depth, particle.lat, particle.lon]

class CustomParticle(JITParticle):
    on_land = Variable('on_land', to_write=True)
    age = Variable('age', dtype=np.float32, initial=0, to_write=True)
    distance = Variable('distance', dtype=np.float32, initial=0)
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False, initial=attrgetter('lon'))
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False, initial=attrgetter('lat'))

this runs fine (i do get a lot of debug info errors that causes the model to freeze with a greater number of particles - is there any way to turn these off? i have tried supressing debug output with pset.kernel_options = {'verbose': False}:

Correct cell not found for (82.550205, 5.437869) after 1000000 iterations
Debug info: old particle indices: (yi, xi) 0 253
            new particle indices: (yi, xi) 0 253
            Mesh 2d shape:  424 324
            Relative particle position:  (xsi, eta) 9.2866782888461774e-01 -8.1678565445934798e-02

however, when i look at the outputs something funny is happening. the particles seem to still be moving over land (see image):

and looking at the variable 'on_land' in the output i see some strange numbers - i am expecting a binary, with 1 denoting particles on land but that is not what i'm seeing in the output.

data:

 on_land =
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0.4727165, 0.7655728, 0.7296033, 0, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2497991,
    0.7931532, 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,

I have tried changing my function to:

def Sample_land(particle, fieldset, time):
    particle.on_land = fieldset.landmask.data[time, particle.depth, particle.lat, particle.lon]

but apparently this isn't supported. Please could you provide some insight into what's going on here?

[erikvansebille]

Dear @mollyjames2, two quick comments:

1. The landmark is a tracer field, for which the default interpolation method is 'linear', so can be any value between 0 and 1. If you want the integers only, you can use 'nearest' interpolation, see https://docs.oceanparcels.org/en/latest/examples/tutorial_interpolation.html
2. While you sample the land (particle.on_land) I don't see that you actually do anything with that. So unless you write a Kernel to tell Parcels that particles should not move onto land, they will still be able to go onto land. Parcels itself has no 'concept' of land, in that sense

[JamiePringle]

are you using the data on the native grid? If you are getting it from Copernicus, they can sometimes distribute it on the A-grid, which will screw things up. The easy way to tell is to see if the horizontal locations of the U and V data points are the same. If they are, you are using the interpolated grid.

[mollyjames2]

i just checked and yes they are the same. Could you please advise how to handle this data given that it's interpolated? thank you

[JamiePringle]

There is no real way to "un-interpolate" it. However, if this is Copernicus data, if you poke arround you can get the c-grid data. You might have to bug their help desk. Otherwise, whoever ran it must have the original data. Sorry I have no more helpful advice. Jamie

[mollyjames2]

Hi,
i plotted the grid points reading in assuming an A-grid and C-grid following the grid indexing tutorial at https://docs.oceanparcels.org/en/latest/examples/documentation_indexing.html

#load the modules for PARCELS
%matplotlib inline
from parcels import Field, FieldSet, ParticleSet, Variable,ScipyParticle, JITParticle, AdvectionRK4,AdvectionDiffusionM1, DiffusionUniformKh, StatusCode
import numpy as np
import math
from datetime import datetime, timedelta
from operator import attrgetter
import warnings
import pandas
import itertools
import os
from datetime import timedelta
import matplotlib.pyplot as plt
from glob import glob

warnings.filterwarnings('ignore')

#set up the velocity fields in a FieldSet object
dataset_folder = "/data/proteus1/scratch/moja/projects/SANH/Leak1H"

# 1.
ufiles = sorted(glob(f"{dataset_folder}/SRIL34_1HL_*1h_grid_U.nc"))  # the name of the file containing the zonal velocity
vfiles = sorted(glob(f"{dataset_folder}/SRIL34_1HL_*1h_grid_V.nc"))  # the name of the file containing the meridional velocity
wfiles = sorted(glob(f"{dataset_folder}/SRIL34_1HL_*1h_grid_W.nc"))  # the name of the file containing the vertical velocity
mesh_mask = f"{dataset_folder}/coordinates.nc" # the name of the file containing the grid specifications for the hydrodynamic data       

filenames = {"U": ufiles, "V": vfiles, "W": wfiles}

# 3.
variables = {"U": "uoce", # the name of the zonal velocity within the files SRIL34_1HL_*1h_grid_U.nc
             "V": "voce", # the name of the meridional velocity within the files SRIL34_1HL_*1h_grid_V.nc
             "W": "woce"} # the name of the vertical velocity within the files SRIL34_1HL_*1h_grid_W.nc

dimensions = {
    "U": {
        "lon": "nav_lon",
        "lat": "nav_lat",
        "depth": "depthu",
        "time": "time_counter",
    },
    "V": {
        "lon": "nav_lon",
        "lat": "nav_lat",
        "depth": "depthv",
        "time": "time_counter",
    },
    "W": {
        "lon": "nav_lon",
        "lat": "nav_lat",
        "depth": "depthw",
        "time": "time_counter",
    },
}

fieldsetA = FieldSet.from_netcdf(filenames, variables, dimensions)



filenames = {
    "U": {"lon": mesh_mask, "lat": mesh_mask, "depth": wfiles[0], "data": ufiles},
    "V": {"lon": mesh_mask, "lat": mesh_mask, "depth": wfiles[0], "data": vfiles},
    "W": {"lon": mesh_mask, "lat": mesh_mask, "depth": wfiles[0], "data": wfiles},
}

c_grid_dimensions = {"lon": "glamf", #define the names of lat and lon from the coordinate file (C-grid)
                     "lat": "gphif",
                     "depth": "depthw",
                     "time": "time_counter"} 




dimensions = {"U": c_grid_dimensions,
              "V": c_grid_dimensions,
              "W": c_grid_dimensions}
# 5.
fieldsetC = FieldSet.from_nemo(
    filenames, variables, dimensions, netcdf_decodewarning=False
)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa

fig, ax = plt.subplots()
nind = 3
ax1 = ax.plot(
    fieldsetA.U.grid.lon[:nind, :nind],
    fieldsetA.U.grid.lat[:nind, :nind],
    ".r",
    label="U points assuming A-grid",
)
ax2 = ax.plot(
    fieldsetA.V.grid.lon[:nind, :nind],
    fieldsetA.V.grid.lat[:nind, :nind],
    ".b",
    label="V points assuming A-grid",
)

ax3 = ax.plot(
    fieldsetC.U.grid.lon[:nind, :nind],
    fieldsetC.U.grid.lat[:nind, :nind],
    ".k",
    label="C-grid points",
)

ax.legend(handles=[ax1[0], ax2[0], ax3[0]], loc="center left", bbox_to_anchor=(1, 0.5))
plt.show()

as shown it appears that my grid is definitely interpolated. with this in mind, would it be appropriate to read in the read with the FieldSet.from_netcdf() function as opposed to Fieldset.from_nemo()

[JamiePringle]

Yes. You will still have issues with with the model not respecting the land boundary, but I bet it would improve things greatly to use from_netcdf becasue then the u and v velocities will be in physically and numerically consistent locations. It is certainly, I think, more physically correct.

Jamie