opendrift simulation with high resolution coastlines by using local nc file

[Morhaf-ALjber]

Dear all, I'm using open drift to run a simulation for a particle tracking in a bay using the water currents resulted from Schism model. I have the binary map of the bay and I want to used it in the simulation because the resolution is higher than those provided by the code. I imported the nc file map which contains longitude, latitude, and land_binary_mask using, reader_landmask = reader_netCDF_CF_generic.Reader('land_mask.nc'), to my code and I changed the debugging source in the line number 83 in "reader_global_landmask.py" to the directory that contains the file that I want to use. Still I'm receiving the resulted plotted map from GSHHG.
How I can obtain my resulted using the high resolution map I have.

[gauteh]

Hi, do you have the chance to share the nc file with the landmask or do a ncdump -h of the landmask file? Please also submit the logs and the essential lines you use for setting up the simulation.

[Morhaf-ALjber]

Dear Gauteh, thank you for your reply.
Apparently, I made a mistake in plotting so the land_binary_mask that I inserted was not used in the plot, now it is working.
The photo attached here is the domain of the simulation identified by the land mask, which is the same domain used in the schism simulation. I also included here the code that I used.

reader_landmask = reader_netCDF_CF_generic.Reader('~/Desktop/opendrift/land_mask.nc')
o.add_reader(reader_landmask, variables='land_binary_mask')
from matplotlib.colors import ListedColormap
import cartopy.feature as cfeature
cmap = ListedColormap((cfeature.COLORS['land'],
cfeature.COLORS['water']))
print(o)
o.plot(background ='land_binary_mask', hide_landmask=True, cmap=cmap, filename = 'test_3.png')

I wonder if I can seed the elements in the river upstream that is shown in the north part of the bay. I do not know if opendrift will use my input nc file as the domain or it will keep the source from GSHHG, since the GSHHG does not include the river and some other water passages. So will opendrift will simulate the particles in the river upstream ?.

Thank you

(https://www.dropbox.com/s/r2f6xh5g2gvb3in/schism_.py?dl=0)
https://www.dropbox.com/s/x6cj4x66oucragu/land_mask.nc?dl=0

[knutfrode]

Hi,
OpenDrift will use your custom landmask for the simulations if you add it as above, and additionally disables the automatic GSHHG landmask
o.set_config('general:use_auto_landmask', False)

Se a demonstration of this here:
https://opendrift.github.io/gallery/example_model_landmask.html

[Morhaf-ALjber]

Dear Prof. Knut-frode, Thank you for your reply.
I understood.

[ocgabs]

Hi @knutfrode
From the example_model_landmask I assume that when using the reader_ROMS_native it reads the model landmask, and therefore you just need to disable GSHHG for the simulation to grab the model grid, correct?
I tried this for my own ROMS domain and it didn't work, even thought mask_rho, mask_v, mask_u are there in the ROMS netcdf. I am thinking to add the landmask explicitly like @Morhaf-ALjber describes... which ROMS mask variable should I use?
Thanks in advance for your help!

[knutfrode]

@ocgabs
In your case OpenDrift should be able to extract the landmask from the ROMS file.
Is land_binary_mask included when you print your reader?
Can you eventually make available a small sample file, so that I can have a look?

[ocgabs]

Thanks for your reply @knutfrode !
When I print the reader the land mask seems to be there, actually it appears twice (see below), however when a run finishes I get:
Model: OceanDrift (OpenDrift version 1.11.1)
0 active Lagrangian3DArray particles (0 deactivated, 0 scheduled)

---

Environment variables:

Readers not added for the following variables:
land_binary_mask
ocean_mixed_layer_thickness
ocean_vertical_diffusivity
etc...

---

Reader: roms native
Projection:
None
Coverage: [pixels]
xmin: 0.000000 xmax: 271.000000 step: 1 numx: 272
ymin: 0.000000 ymax: 316.000000 step: 1 numy: 317
Corners (lon, lat):
(152.85, -25.12) (162.22, -28.16)
(147.10, -38.50) (158.03, -41.55)
Vertical levels [sigma]:
[-0.97146552 -0.91560345 -0.86135057 -0.8087069 -0.75767241 -0.70824713
-0.66043103 -0.61422414 -0.56962644 -0.52663793 -0.48525862 -0.44548851
-0.40732759 -0.37077586 -0.33583333 -0.3025 -0.27077586 -0.24066092
-0.21215517 -0.18525862 -0.15997126 -0.1362931 -0.11422414 -0.09376437
-0.07491379 -0.05767241 -0.04204023 -0.02801724 -0.01560345 -0.00479885]
Available time range:
start: 2009-12-02 00:00:00 end: 2009-12-06 21:00:00 step: 3:00:00
40 times (0 missing)
Variables:
ocean_vertical_diffusivity
sea_floor_depth_below_sea_level
latitude
latitude
latitude
latitude
longitude
longitude
longitude
longitude
land_binary_mask
land_binary_mask
ocean_s_coordinate_g2
ocean_s_coordinate_g2
sea_water_salinity
surface_downward_x_stress
surface_downward_y_stress
sea_water_temperature
x_sea_water_velocity
y_sea_water_velocity
upward_sea_water_velocity
sea_surface_height
sea_water_speed - derived from ['x_sea_water_velocity', 'y_sea_water_velocity']

Hopefully you can download a sample file here: https://www.dropbox.com/scl/fi/1py75fvsa068e34kenwnm/part_ocean_fwd001_blend_04883.nc?rlkey=jx6mo80tx4fwhfczp43249564&dl=0

Thanks a lot for looking into this. It means the world!

[ocgabs]

Hi again @knutfrode,
I hope you can help me because it seems that I have serious issues... perhaps they are all related to the land mask not being properly read... It seems that loads of particles leave the domain and keep going! In this plot I overlayed the model bathymetry to mark the extend of the ROMS model on top of the end position of particles in magenta stars.
One particle that I removed, to help the visualization jumped over the globe and ended up at longitude -150!
I am including my code below. Am I doing something wrong? Even if the time step is too long, I would expect just to jump out of the domain and then stop, and not getting this far :/
I will go back to advection only experiments and see if its better...

#!/usr/bin/env python
"""
ROMS native reader

"""

import numpy as np
from datetime import datetime
from datetime import timedelta
import matplotlib.pyplot as plt
from glob import glob

from opendrift.readers import reader_ROMS_native
from opendrift.models.oceandrift import OceanDrift
from opendrift.readers import reader_global_landmask, reader_netCDF_CF_generic #, reader_netCDF_CF_unstructured # for corners

data_path = '/scratch/public/ca5240/2018/'

Reader Ocean Current

reader_ROMS = reader_ROMS_native.Reader(sorted(glob(data_path+'/ocean_fwd001_blend_0*.nc')))

print(reader_ROMS)
reader_ROMS.Dataset

print(reader_ROMS.start_time)

print(reader_ROMS.end_time)

corners = reader_ROMS.xy2lonlat([reader_ROMS.xmin, reader_ROMS.xmin, reader_ROMS.xmax, reader_ROMS.xmax],
[reader_ROMS.ymin, reader_ROMS.ymax, reader_ROMS.ymax, reader_ROMS.ymin])
corners

Set particle simulation

o = OceanDrift(loglevel=20) # Set loglevel to 0 for debug information
o.add_reader([reader_ROMS])

o.set_config('drift:horizontal_diffusivity', 1000)
o.seed_elements(lon=154, lat=-26, radius=5000, number=40000, time=reader_ROMS.start_time)
o.seed_elements(lon=150.963, lat=-35.403, radius=5000, number=40000, time=reader_ROMS.start_time)
duration=timedelta(days=365)

RUN!

o.run(time_step=1800,duration=duration, time_step_output=24*3600, outfile='2sHD5000_ODtestNCIv2y2018.nc') `

[knutfrode]

You have quite a bit of diffusivity (1000 m2/s), so I guess particles leave the domain by this mechanism.
For the OceanDrift models wind and current are not strictly required, allowing to to sensitivity tests. If you want particles to be deactivated whenever wind and/or currents are missing (as e.g. for OpenOil and Leeway), you can specify:
>>> o.set_config('environment:fallback:x_wind', None)

I also recommend setting loglevel to 0. Then you get a detalied log, that will tell you everything that happens during simulation.
E.g. it will tell why seafloor depth was not used.

[ocgabs]

Thanks a lot for your answer @knutfrode
That works so particles don't leave the model domain.
Also, I have now managed to turn off the automatic land-mask, perhaps the problem was just the order of the config commands.

I can't figure the syntax to address this warning: Data block from roms native not large enough to cover element positions within timestep. Buffer size (5) must be increased. See Variables.set_buffer_size.
Can you please share an example?

Also, it seems that in previous versions the run command would take the option corners, but not anymore? How can I include corners? it seems that they are needed for some of the plot functions to work...

Thanks a lot for all your help and for making this amazing code. I am now exploring the plotting and analysis tools. It will be so great if I can reproduce many of the figures on the gallery!

[knutfrode]

The warning about buffer_size may be a bit misleading as it also happens when particles are close to the reader domain, so that it is not possible to read a block large enough to cover future propagation. This warning will be removed in the near future, but for now you can neglect it.

I believe corners was never an option to the runmethod, just for the plotting, as in e.g. this example:
https://opendrift.github.io/gallery/example_grid.html

[ocgabs]

Thanks a lot for your response! My concern about the buffer_size warning is because I get the following when running:
11:17:53 WARNING opendrift.readers.basereader.structured:286: Data block from roms native not large enough to cover element positions within timestep. Buffer size (6) must be increased. See Variables.set_buffer_size.
11:17:53 WARNING opendrift:133: Still NaN-values after 10 iterations, exiting!
11:17:53 WARNING opendrift:133: Still NaN-values after 10 iterations, exiting!
11:17:54 WARNING opendrift.models.basemodel:2276: Missing variables: ['x_sea_water_velocity']
But then it keeps running... so I am not sure if there is a problem...