eliocamp / metR

Tools for Easier Analysis of Meteorological Fields
https://eliocamp.github.io/metR/
139 stars 22 forks source link

streamlines / Trajectory example #149

Open danstrobridge-Weston opened 3 years ago

danstrobridge-Weston commented 3 years ago

Hi there, I'm interested in plotting trajectories of several particles based on a steady-state groundwater flow field in a small area that is projected into state plane coordinates. I suspect that this would be a very simple use case of your Trajectory() function. I'm hoping you might provide an example that I could use since I don't see anything in the help file or on the git page. I've got a groundwater elevation raster that I've converted to data.table (DT) format. I don't have the physics background of a metereologist but can easily calculate slope and/or aspect of the raster cells and add those to the DT, if necessary. In your formula notation, I gather dx and dy are the derivatives with respect to x and y and the x and y are the coordinates of the raster cells themselves? If true, then I can easily add dx and dy to the DT that I'm assuming would be employed via the "data" option. I see that the function works on just one particle, so I'd need to purrr::map on a data frame of desired particle starting points... unless you've already got a wrapper utility that does this.

At first, I tried following the example in the help of geom_streamline() but then saw dx and dy were calculated with GeostrophicWind, which is not applicable for my use. I would indeed love to plot the streamlines via ggplot on top my rasterVis based plot object. However, I do also want to be able to export the trajectory to a line shapefile for use in GIS as well.

eliocamp commented 3 years ago

If what you want to compute are trajectories of a static velocity field, then you can use geom_streamline(). Trajectory() is for time-varying velocity fields.

Assuming that velocity is directly proportional to negative slope, then this would be an example:

library(metR)
library(ggplot2)
library(data.table)
data(volcano)
volcano <- as.data.table(melt(volcano))

volcano[, c("dx", "dy") := Derivate(value ~ Var1 + Var2)]

ggplot(volcano, aes(Var1, Var2)) +
    geom_contour_fill(aes(z = value)) +
    geom_streamline(aes(dx = -dx, dy = -dy), n = 20, res = 2) +
    coord_equal()

Created on 2021-08-06 by the reprex package (v2.0.0)

The first few lines are to compute the necessary data, you don't need to use data.table if you don't want you. What you need to have is a table with x and y locations with their respective dx and dy values .

fleetboat commented 2 years ago

I have a time-varying velocity field (u and v velocity components as functions of position (lat, long) and time (t) and would like to use the Trajectory function to determine the path of a fluid parcel originating at a specific point in the field. It would seem that this is exactly the purpose of the Trajectory function in metR, but the documentation neither explains what is returned by the function nor provides an example of its use. Do I understand correctly that for my data, dx and dy would be the derivatives (finite differences) of the velocity components with respect to time. The input data frame then would have variables xposition, yposition, u(at xposition, yposition), v(at xposition, yposition), du (at xposition, yposition), dy(at xposition, yposition), and t so would be a frame with dimension [gridXsize gridYsize t rows, and 7 columns]?

Thank you

eliocamp commented 2 years ago

Sorry for the lack of documentation. The input field would be the velocity components at each grid point and each time step. You don't need the acceleration (derivatives of the velocity components).

The result is a data.frame with the x and y position of each starting particle at each time along with its velocity.

This would be an example (although the resulting trajectory is rather boring because I could not come up with an interesting field.

library(ggplot2)

x <- -10:10
y <- -10:10
times <- 1:10

grid <- expand.grid(x = x, y = y, time = times)

grid$dx <- (-(times-5)^2 + 2.5^2)/10
grid$dy <- (-(times-5)^2 + 2.5^2)/10

trajectory <- metR::Trajectory(dx + dy ~ x + y + time, x0 = c(5, 3), y0 = c(2, 1), data = grid)

ggplot(trajectory, aes(x, y)) +
    geom_path(aes(group = id))

Created on 2021-09-04 by the reprex package (v2.0.0)

fleetboat commented 2 years ago

Thank you. That is very helpful. A few notes that may help other users: First, the dx, dy units have to be compatible with the x, y units. Without thinking about it, I used geographic coordinates for x and y, and m/s for dx, dy in my first test. That failed. When I converted the grid to distance units (meters) the trajectories were correct. Second, the input data frame can get very large, very quickly especially when geographic positions (as opposed to grid indices) are used for the grid. I was testing with a 61x30 point grid that had 1400 time steps. My R installation choked on that.