OpenDrift / opendrift

Open source framework for ocean trajectory modelling
https://opendrift.github.io
GNU General Public License v2.0
249 stars 120 forks source link

seed_elements issue #478

Closed trondkr closed 3 years ago

trondkr commented 3 years ago

Hi,

I wonder if I am doing something wrong here or if there is a non-logical issue with the seed_elements function in basemodel.py. If you have a list of longitudes and latitudes and you want to seed 100 particles in a radius around each of these locations between start _date and end_date I would assume that you would write this:

o.seed_elements(lon=[lon_list],
                        lat=[lat_list],
                        number=100,
                        radius=radius,
                        time=[start_date, end_date]) 

However, OpenDrift complains about length of lons not being equal to number (raise ValueError('Lon and lat have length %s, but number is %s' % (len(lon), number))). Looking into basemodel.py I see that the code goes through these lines:

        if len(lon) > 1:
            if number is not None and number != len(lon):
                raise ValueError('Lon and lat have length %s, but number is %s' % (len(lon), number))
            number = len(lon)

Is this correct behavior? I see no reason why the particle number would have to be equal to the length of station longitudes? I can see that the list of numbers (if provided) would have to be equal but if there is a constant number of particles to be released at each location this seems wrong.

I am using Opendrift 1.4.2 installed using pip.

Curious what you think.

Thanks, Trond

knutfrode commented 3 years ago

Hi Trond,

I wonder if your expected behaviour was actually happening until a few months ago. But the old seed_elements was very slow due to recursivity and complexity, and there was also some ambiguity: E.g. if your lon/lat_list has 10 points, and you ask for 100 elements, it could mean either 100 elements at each point (as you say you expect), but it could also mean 10 elements for each of the 10 positions.

So now you would have to use either a loop (which is less of a problem since each seed-call us much faster than before), or you can prepare your arrays to have equal length by e.g. meshgrid. A bit more cumbersome sometimes, but you have more control without ambiguity.

Note that time may also have length equal to the number of elements (or seed locations).

trondkr commented 3 years ago

Hi and thanks for the explanation. I was probably expecting the old behavior but I understand the reasoning for changing this. Does this also mean that if I seed 1000 particles at one location with start and stop dates (in a list), the particles will be distributed in time and seed a total of 1000 particles, or will 1000 particles be seeded at each timestep?

Thanks for helping. Enjoy the holidays :) T.

knutfrode commented 3 years ago

Yes, only time is still treated the way that if it is a list of two elements, the total number of elements (e.g. 1000) will be seeded linearly/evenly within this time interval. All other numbers (lon, lat, radius, model-specific-seed variables...) must be either a scalar (same for all elements) or an array of length equal to the number of elements. If number is not given explicitly, it is either automatically the length of any arrays (e.g. lon/lat or other variables), or it is taken from the config setting seed:number

I believe the updated docstring should be correct, though admittably a bit hard to comprehend: https://opendrift.github.io/autoapi/opendrift/models/basemodel/index.html#opendrift.models.basemodel.OpenDriftSimulation.seed_elements

trondkr commented 3 years ago

Thank you!

CaroMedel commented 2 years ago

Hi,

I have the same problem, I want to seed particles in different points. It works well if I don't define number, but I want to define the number. I tried this and it didn't work:

num_steps = 73 time_step = timedelta (hours = 1 )

for i in range (num_steps ): o .seed_elements (lon=[-75.0531, -75, -75.4], lat=[-43.2, -43.7, -42.5], number = [10,10,10], time =time +(i *time_step ))

Could you help me? Thanks

knutfrode commented 2 years ago

Hi,

number must be an integer, and not an array. But you may loop over the positions and seed over a time interval instead of looping over time steps:

from datetime import datetime, timedelta
from opendrift.models.oceandrift import OceanDrift

o = OceanDrift()

lons = [-75.0531, -75, -75.4]
lats = [-43.2, -43.7, -42.5]
start_time = datetime(2022, 2, 1)
end_time = datetime(2022, 2, 2)

for lon,lat in zip(lons, lats):
    o.seed_elements(lon, lat, number=100, time=[start_time, end_time])

o.set_config('environment:constant:x_sea_water_velocity', 1)
o.set_config('environment:constant:y_sea_water_velocity', 1)
o.run(duration=timedelta(days=3))

o.animation(fast=True)

This is slightly different in that elements are spread uniformly over the time interval instead of in groups of 10 every hour. This might however be better(?)

If you prefer to seed them in groups, you can prepare arrays with meshgrid so that number is an integer for each seeding call.