"Missing variables: ['land_binary_mask']" at the end of simulation?

[IreneNA95]

Hi everyone!

I am performing particle advection experiments using the Opendrift code and the SHYFEM native readers. Following indications from https://opendrift.github.io/gallery/example_shape_landmask.html#use-a-shapefile-as-landmask, I am using a shapefile as customized landmask. The code that I insert in my simulation is:

from opendrift.readers import reader_shape
....
reader_coast = reader_shape.Reader.from_shpfiles('coast.shp')
o.add_reader(reader_coast)
o.set_config('general:use_auto_landmask', False)  # Disabling the automatic GSHHG landmask

Actually, when I run the simulation, it is confirmed that the imposed landmask is properly recognized: 00:31:00 INFO opendrift.models.basemodel:1689: Using existing reader for land_binary_mask However, after few time steps, I obtain the following warning:

00:38:03 INFO    opendrift.models.basemodel:2870: 2018-07-01 15:00:00 - step 31 of 2880 - 1200 active elements (0 deactivated)
00:38:15 WARNING opendrift.models.basemodel:3133: Missing variables: ['land_binary_mask']
00:38:16 INFO    opendrift.models.basemodel:2870: 2018-07-01 15:30:00 - step 32 of 2880 - 1199 active elements (1 deactivated)

And then, at the end of simulation, a large quantity of the same warning "Missing variables: ['land_binary_mask']"is displayed.

When I plot the trajectories with the coastline as background, I see that some of my particles overpass the coastline (even overpassing the available hydrodynamic grid).

Could someone give me some advice on how to fix this issue?

[knutfrode]

Hi,

It seems that there is only one particle that has missing variable land_binary_mask for this particular timestep. So it might seem that some particles are leaving the domain of the shape-reader?

Note that if calculation time step is too large (depending on application), particles may "jump" over land pixels. So if this is small and localized area, you'd might like to reduce calculation time step to e.g. 5 minutes to avoid this problem. But this should not be related to the "missing variables" message.

When plotting, you should see the discarded particles in gray. Does this seem to happen at the boundary of the landmask reader? Can you post such a plot here? And maybe you can also post your full invoking script?

[IreneNA95]

Hi, Thank you really much for your fast reply. Yes, actually my time step is 30 minutes. I will try to decrease it and see if the problem persists.

I post a plot of trajectories here:

And I post my full code here:

#!/usr/bin/env python
"""
SHYFEM: Using model input from unstructured grid
================================================
"""
#---- LOADING PYTHON LIBRARIES ----#

# Importing useful libraries of python:
import numpy as np
from datetime import timedelta
import matplotlib.pyplot as plt
import netCDF4 as nc
import scipy.io
# Importing opendrift readers and model:
from opendrift.readers.unstructured import shyfem      # Readers from SHYFEM.
from opendrift.readers import reader_netCDF_CF_generic # Other useful and recommended readers for variables not readed with shyfem functions.
from opendrift.readers import reader_shape
from opendrift.models.oceandrift import OceanDrift

# Creating a simulation object by calling the imported model, with minimum level of verbose (20)
o = OceanDrift(loglevel=20)

#---- LOADING DATA ----#

# Using shyfem readers to read shyfem model outputs (=inputs of lagrangian algorithm).
shyfem = shyfem.Reader(filename='inputopendrift/JAS/shyfem_2018*.nc')
o.add_reader(shyfem)
print(shyfem) # Print the read variables from SHYFEM file

# Loading coast
o.set_config('general:use_auto_landmask', False)  # Disabling the automatic GSHHG landmask
reader_coast = reader_shape.Reader.from_shpfiles('coast.shp')
o.add_reader(reader_coast)

# Loading release locations
mat = scipy.io.loadmat('INPUT/valid_release_adriatic_res1.00.mat')
cx = mat['CX'] # center longitudes of cells
cy = mat['CY'] # center latitudes  of cells

#---- DEFINING SEEDING PARAMETERS ----#

# Defining depth (in meters)
z = -1

# Defining radius (in meters)
r = 12000

# Defining number of particles
g = 50
n = np.size(cx)*g # 50 particles per cell

x = np.reshape(np.tile(cx,(1,g)),(-1,1))
y = np.reshape(np.tile(cy,(1,g)),(-1,1))

print(np.size(x))
print(np.size(y))
print(n)
# Seeding n elements at defined lon,lat (x,y) positions, depth (z) and time.
o.seed_elements(lon=x, lat=y, radius=r, number=n,z=z, time=shyfem.start_time)

# Interaction with the coastline. 
o.set_config('general:coastline_action', 'previous') 

# Lagrangian propagation scheme used for the run. Runge-kutta-4th order: most used method for lagrangian dispersion.
o.set_config('drift:advection_scheme', 'runge-kutta4') 

o.set_config('vertical_mixing:TSprofiles', True)

# Name of output file
fo = 'LPT_Adriatic_30mnJAS.nc' # LPT: Lagrangian Particle Tracking

#After initialisation and configuration, a model run can be started by calling the function o.run():
o.run(time_step=1800, duration=timedelta(days=60), outfile=fo)

#---- PRINT AND PLOT RESULTS ----#
print(o)
#o.plot(fast=True,linecolor='z') # this will color lines according to depth (z):
#o.animation(color='z', markersize=5, corners=[12.2, 12.6, 45.1, 45.5],filename='animation.mp4')

Thanks again,

[knutfrode]

Hi,

This figure shows that the discarded elements are found only at the southern part, which is probably outside of your shapefile, and thus it is then expected that they should be deactivated, as coastline is then unknown.

It should be noted that the landmask shown on the figure is not your own provided (which is used during simulation), but the GSHHG landmask. There seem to be a small shift among these, as trajectories are generally confined to the sea, but with some examples of onland propagation in eastern shore. The GSHHS geolocation is pretty accurate; could it be that the one you are using is slightly offset somehow?

Btw, if you add lscale='f' (or 'h') to plot and animation, you will use the GSHHG version with full (or high) resolution. I.e. with more precision than what is used for your plot. But plotting then also takes more time.

[IreneNA95]

Hi, I have plotted (in MATLAB) the GSHHS in the full mode (with m_gshhs_f function) and the coastline that I am providing to the simulation, and it seems not displaced, so it might be due to a different reason.

Anyway, is there a way in which I can use the GSHHS in the "full resolution" mode as a landmask in my simulation?

Thank you again for your availability.

[knutfrode]

For the simulation, the full resolution GSHHG landmask is always used. It is only for the plotting that a lower resolution might be used to make it faster. For a large area (like this one), 'h' or 'i' might have more than enough details for the plot. But if not chosen, plot/animation methods will chose an automatic precision level.

Thus you can probably skip adding the customized landmask for this area (and instead use the built in GSHHG mask).

[IreneNA95]

Hi, Running the simulation with the automatic GSHHG mask (maintaining the time step of 30 minutes) resulted in a much faster simulation (besides the warning at the beggining of simulation: Adding a dynamical landmask with max. priority based on assumed maximum speed of 3 m/s. Adding a customised landmask may be faster...) and, as expected, solved the problem of the warning Missing variables: ['land_binary_mask'] during and at the end of simulation.

In any case, I will follow your suggestion of decreasing the time step to 5 minutes. Thank you really much for your fast replies and good advices.

[knutfrode]

Very good. I guess time step of 15, 30 or even 60 minutes should probably be fine at this scale.

Then closing this issue.