Residence time?

[AndresSepulveda]

I would like to use OpenDrift to study the residence time of a semi-enclosed area. I would need to

1. Prevent particles from stranding
2. Fill the area with as many particles as possible
3. Discard those leaving a defined boundary

Any hinds on how (if possible) to define this here? Thanks

[knutfrode]

Hi,

This should be possible to achieve somehow.

1. Presently, all OpenDrift models deactivate particles hitting the coastline (except for WindBlow).
   I have plans to later make other options available, e.g. to keep particles at the coast until they are moved offshore at a later time.
   In the meantime, the best is perhaps to try to avoid letting particles hit the coast, by using as small (calculating) time_step as possible (time_step_output may still be large). Assumedly, your ocean model will have vanishing, or at least small, onshore current component.
   It may also help to use Runge-Kutta scheme instead of the default Euler scheme:
   o.config['drift']['scheme'] = 'runge-kutta'
   You may also try to run without a landmask, but I think particles will then simply be stuck anyway when leaving the "water pixels" in your ocean model:
   o.fallback_values['land_binary_mask'] = 0

2. To fill the area with particles, you may use the seeding method seed_within_polygon:
   https://github.com/knutfrode/opendrift/blob/master/opendrift/models/basemodel.py#L897
   Input is arrays of lon and lat of a polygon (your semi-enclosed area), and the number of particles.

3. One way to discard particles leaving a defined boundary, is to let your Basemap reader not cover this area. Note that it is possible to define Basemaps using any map projections, including rotated. In this case fallback_values for land_binary_mask must not be set to zero as above, because particles will then not be deactivated when leaving the "map".
   Alternatively, you may make your own modified version of the OceanDrift module, which would give you total freedom:
   https://github.com/knutfrode/opendrift/blob/master/opendrift/models/oceandrift.py
   E.g. to remove all particles moving north of, say -36N, you may add the line:
   self.deactivate_elements(self.elements.lat > -36, reason='leaving')

[AndresSepulveda]

Great, that solves 2) and 3). As for 1), there could be also two other options

A) For the particles to be able to bounce on the coast.
B) To set to zero the velocity component that is perpendicular to the coast.

In both cases one would need to know it the nearby cells are land. Is there an option to obtain that information for each particle, at each time step?

[knutfrode]

Yes, those are two other relevant options.

Regarding "land in nearby cells":
The OpenDrift core (and thus the various models/subclasses) operates in a spatial reference system ("map projection") where each particle has a position (x,y), corresponding to a given (lon,lat). But there is no discretisation in space here, so no delta_x and delta_y, and no "nearby cell".

The readers, however, are responsible to obtain forcing data (current, temperatures, wind etc) from various external sources, which is most often output from an ocean model etc, with data on a discrete, regular grid. The readers then also interpolate this data onto the particle positions, and provide this data to the actual model, which knows nothing about how the data was obtained.
The "Basemap reader", which normally provides the landmask towards which particles are stranding, is also a regular reader, which can be replaced with any other reader which provides "land_binary_mask" (e.g. often provided by ocean models). The Basemap reader does however not obtain data from a discretised grid, but returns 1 for each particle/position which is inside any land contours from the global GSHHS dataset, and 0 otherwise (ocean). This is simple and fast, but there is a danger that particles may "jump" over islands if the calculation time step is too large.
A more sophisticated interaction with land could be implemented in various ways, but it is probably a good ide to build something around the Basemap reader, as this is a global, high resolution dataset, which is always available. Thus when a particle is moved from A to B (due to currents, wind etc), a mechanism is needed to find the land-crossing point, and then calculate a new position by either sticking or deactivating the particle here, or bouncing or sliding, as you suggested. Such calculation would probably take significant calculation time, as it must be done/checked for all particles at all time steps. It should probably be parallelised code to be manageable.
Unfortunately, I will not be able to prioritise implementing this in the near future. (For my own main applications, oil and search&rescue, stranding towards the coast is desirable!) But if anyone would be interested to look at implementing non-stranding options in OpenDrift, I would be happy to assist.

[AndresSepulveda]

Thaks for the reply. I'll try to find someone better at coding than me to work on those ideas.

I understand from your comment that the "land_binary_mask" could be read from the forcing file, a ROMS file in my case. If so, do you have an example of how this is done?

[knutfrode]

You can check from commandline $ readerinfo.py <ROMS_file> or from Python >>> print <ROMS_reader> if the forcing file contains the variable "land_binary_mask".
In that case, if you add this reader and no basemap-reader (or if you add the ROMS-reader before the basemap-reader), the landmask from the ROMS file will be used during the run.

[AndresSepulveda]

Strange.. the file has a variable with that description

    float mask_rho(eta_rho, xi_rho) ;
            mask_rho:long_name = "mask on RHO-points" ;
            mask_rho:option_0 = "land" ;
            mask_rho:option_1 = "water" ;
            mask_rho:standard_name = "land_binary_mask" ;
            mask_rho:coordinates = "lat_rho lon_rho" ;

but the script readerinfo doesn't recognizes it

andres@loki:~/opendrift/tests/test_data/Piti$ ../../../opendrift/scripts/readerinfo.py roms_his_v6T.nc
Testing /home/andres/opendrift/opendrift/readers/reader_netCDF_CF_generic.pyc...
...not applicable.
Testing /home/andres/opendrift/opendrift/readers/reader_ROMS_native.pyc...

Reader: roms native
Projection:
None
Coverage: [pixels]
xmin: 0.000000 xmax: 57.000000 step: 1 numx: 58
ymin: 0.000000 ymax: 234.000000 step: 1 numy: 235
Corners (lon, lat):
(-72.84, -43.65) (-72.80, -43.65)
(-72.84, -43.77) (-72.80, -43.77)
Vertical levels [sigma]:
[-0.94999999 -0.85000002 -0.75 -0.64999998 -0.55000001 -0.44999999
-0.34999999 -0.25 -0.15000001 -0.05 ]
Available time range:
start: 2000-01-31 00:00:00 end: 2000-03-01 00:00:00 step: 1:00:00
721 times (0 missing)
Variables:
sea_floor_depth_below_sea_level
sea_surface_height
x_sea_water_velocity
y_sea_water_velocity
sea_water_temperature
sea_water_salinity

I use ROMS_AGRIF, just in case. Any clue why is not identified?

[knutfrode]

It seems that your ROMS file has the landmask defined on rho-points ("mask_rho"), whereas the ROMS files we have been using (not AGRIF) have landmask defined on psi-points ("mask_psi").

I have now updated the ROMS reader to map both of these onto the CF standard name "land_binary_mask".
Can you please update your OpenDrift version and try again?

[AndresSepulveda]

ok, this works fine. My script is like this now

o = OceanDrift(loglevel=0) # Set loglevel to 0 for debug information

nordic_native = reader_ROMS_native.Reader(o.test_data_folder() +
'Piti/roms_his_v6T.nc')
o.add_reader([nordic_native])

##Landmask (Basemap)
##reader_basemap = reader_basemap_landmask.Reader(
##llcrnrlon=-72.86, llcrnrlat=-43.8,
##urcrnrlon=-72.80, urcrnrlat=-43.64,

resolution='f', projection='merc')
##o.add_reader([reader_basemap, nordic_native])

However it plots an area much larger than the domain, with a coastline (from I-don't-know-where). How can I restrict the area to be plotted within the reader_ROMS_native?

[knutfrode]

Your simulations above are using landmask from ROMS to checking stranding during the run, but for the plotting, Bamemap coast contours are always used.
You may zoom into the figure with the mouse, but you may also add a Basemap reader with customised borders (before plotting) to get plot of your area of interest. This Basemap reader should be added after the ROMS-reader, so that ROMS landmask, and not Basemap landmask, is used for the stranding.

However, even when using ROMS landmask, particles will hit land and become deactivated.
But a very simple option to avoid this, is simply to move any particles hitting land back to their previous positions. A quick example of this is attached, but not tested. With this example, particles should be kept close to shore as long as there is an onshore current component, but be moved offshore as soon as there is an offshore current component at that position. If this works well for you, I will consider adding this as a general alternative to stranding/deactivation in OpenDrift.
no_stranding.zip

[knutfrode]

A very basic non-stranding algorithm is now implemented in OpenDrift, as described in #32
So by specifying o.set_config('general:coastline_action', 'previous'), particles hitting land will simply be moved to their previous offshore position.

This is not well tested yet, but please give it a try.

[knutfrode]

For the ROMS reader there were some problems with the anti-stranding mechanism due to staggere/shifted grids of u and v. This is now improved, although some particles may still be trapped at the coastline.
The different stranding options are now illustrated here.