OpenDrift / opendrift

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

Saving only start and final lat/lon for each particle #582

Closed calquigs closed 3 years ago

calquigs commented 3 years ago

Hi all,

For many of my simulations, the only real pieces I'm interested in are the seed+destination coordinates, and I'm trying to think of a way to get those data and exclude as much extraneous information as possible. Is there already a setting in opendrift to extract those two locations, or is it possible to write a module that saves that data in a new variable?

Thank you, Calvin

gauteh commented 3 years ago

You should be able to slice the array with particle positions to get what you want, is there any reason that doesn't work?

gauteh commented 3 years ago

I think this should give you what you need:

# Running model
lw.run(duration=timedelta(hours=48), time_step=900, time_step_output=3600)

lons, lats = lw.get_lonlats()
print(lons.shape)

start_lon, start_lat = lons[:, 0], lats[:, 0]
end_lon, end_lat = lons[:, -1], lats[:, -1]
calquigs commented 3 years ago

That's very helpful, thank you! My particles aren't all seeded at the same time, so [:,0] and [:,-1] return masked values for most particles, but I can use .compressed() to get unmasked values only for each particle in lons. I can then use this to write to file only what I need at the end of the run rather than write the entire trajectories inside with run(outfile).

This approach still requires that my time_step_output be small enough to capture the start and final location for each particle. Would it be possible to save start/final coordinates somewhere else, and then set the external timestep to something much larger? Not sure how much time this would save, it may not be worth the trouble— just trying to think of more ways to optimize my runs!

knutfrode commented 3 years ago

Hi Calvin,

There is an internal, undocumented method index_of_activation_and_deactivation which can be used like this:

index_of_first, index_of_last = o.index_of_activation_and_deactivation()
lons = o.get_property('lon')[0]
lats = o.get_property('lat')[0]
lons_start = lons[index_of_first, range(lons.shape[1])]
lons_end = lons[index_of_last, range(lons.shape[1])]
lats_start = lats[index_of_first, range(lons.shape[1])]
lats_end = lats[index_of_last, range(lons.shape[1])]

There is however a small issue with seeding location, as you mentioned: with continuous seeding in time, some particles will be seeded/released between two time steps, and the position of these will not be recorded until the end of that (output) time step, when they may have moved slightly from the actual seeding location. For the final position, this is not an issue, as the actual position of deactivation (or the final drifting position) is returned by the method above.

So, as you mentioned, you would need a very small time_step_output to minimize the error of seeding locations. However, you may also recover the exact seeding positions after seeding but before starting the simulation (run) with:

lons_start = o.elements_scheduled.lon
lats_start = o.elements_scheduled.lat
calquigs commented 3 years ago

This has worked very well for me, thank you!

AndresSepulveda commented 3 years ago

Hi Calvin Quigley,

Can you share the final implementation that you used? I am also interested in this kind of features.

Regards,

Andrés

calquigs commented 3 years ago

@AndresSepulveda Hi yes, sorry about that!

Per knutfrode's suggestion, I'm running:

lons_start = o.elements_scheduled.lon
lats_start = o.elements_scheduled.lat

o.run()

index_of_first, index_of_last = o.index_of_activation_and_deactivation()
lons = o.get_property('lon')[0]
lats = o.get_property('lat')[0]
status = o.get_property('status')[0]
lons_end = lons[index_of_last, range(lons.shape[1])]
lats_end = lats[index_of_last, range(lons.shape[1])]
status_end = status[index_of_last, range(lons.shape[1])]

And then just writing to a .txt file with a line containing "lon_start, lat_start, lon_end, lat_end, status_end" for each particle.

AndresSepulveda commented 3 years ago

Just for completeness, this is what I am using

   all = np.column_stack([lons_start, lats_start, lons_end, lats_end, status_end])
   np.savetxt(filename, all, fmt='%.4f')