Friday, 15 July 2016

2D particle tracking in R

How hard is it nowadays to simulate particle movement in the sea and track them in R? Well I discovered that it is not all that difficult. Within a couple of hours I had a simple simulation model running in the North Sea. All you need are time dependent water currents and wind data, both of which can be obtained from marine.copernicus.eu.

The main steps that you will have to follow:

  • load the water currents and wind data
  • specify particle release details
  • specify simulation details
  • loop all simulation time steps and for each step do the following with the particles that are released:
    • for each particle determine the local water current and wind vector at the current time index
    • transform the particle coordinates to a projection in meters
    • move the particles with the water current (m/s) multiplied by the time step duration (s)
    • move the particles with the wind speed (m/s) multiplied by the time step duration (s) and a drift factor (between 0 and 1)
    • move the particles in random directions (to simulate diffusion)
    • transform the particle coordinates back to the original projection
    • save particle information every nth simulation step, where n is any positive integer value

Below you will find the script that does all this (each of the steps above are marked in the script). In this simulation I simulate the release of 5,000 particles with a constant release rate during 5 days and followed them for a couple of days. I found the trickiest bit not the simulation itself, but reading the the NetCDF files such that it could be used in a convenient manner. I create a video of this simulation with ffmpeg:

So what's the practical use of all this? Well obviously you can simulate the trajectory of any floating materials, such as oil. However, you have to keep in mind that the simulation is only as good as the underpinning data and parameterisation. For instance, the water current data used here is not that accurate near estuaries. Of course you can always do multiple runs of your simulation with different parameter settings (or make it even a Monte Carlo simulation). Also, here I only simulate in 2D. But there is nothing that would stop you from modifying it to include the third dimension. What also should be added is some detection that particles have beached. But that should not be very difficult when you have coastline data.

Why in R? Well, commercial particle tracking software is generally expensive. Whereas free software are often rigid and don't allow the user to include custom modifications. In R the simulation is free, and you can take the exact implementation into our own hands... Hopefully, this post will help you get started with your own simulations. Good luck!

Here's the script:

Edit: if your interested, you might want to check out the follow-up post on particle density estimates of this simulation.

5 comments:

  1. Very Helpful post. Thank you very much.
    But I can't find the data used in this code.
    would you let me know the download link?

    ReplyDelete
    Replies
    1. Hi there, thank you for your interest. Both wind and current data were extracted from Copernicus (http://marine.copernicus.eu). Note that the URL occasionally changes, but this is were the current data is now: http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=NORTHWESTSHELF_ANALYSIS_FORECAST_PHYS_004_001_b I think this is where the wind data can be found nowadays: http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=WIND_GLO_WIND_L4_NRT_OBSERVATIONS_012_004

      Note that you do need an account to access the data, but this should be free. Please let me now if you need any more help...

      Delete
    2. Sorry I'm late. I used only the concept you proposed. With my own data, simulation worked pretty well. But.. when particles penetrated in the land, they couldn't get out of the land. If you have an idea about this boundary problem, plz let me know how to deal with it.

      Delete
  2. Hello
    Thank you for this good script to model particle movement in 2D, it will be very helpful once I can get it to work for my data. I am having an issue when I get to this step

    'part.coords <- SpatialPointsDataFrame(part.coords, data.frame(released = rep(F, release.n)))'

    I therefore wanted to check what the 'data.frame' part of this step represents - as I understand it this is to create a dataset with 'FALSE' replicated by the number of particles released. When I try to run this step I get the following error:

    'Error in SpatialPointsDataFrame(part.coords, data.frame(released = rep(F, :
    row.names of data and dimnames of coords do not match'

    Any advice would be greatly appreciated!

    Kind regards

    Jon

    ReplyDelete
    Replies
    1. Update: I managed to fix this issue - it was being caused by some NAs in the part.coords dataframe step, so I just converted these to lat long values and it solved the issue.

      Delete