madpk commented 10 months ago

I run the following code to create a HYSPLIT backward model. However, the _trajectoryplot() function plots only the base map, not the trajectories.

HYSPLIT Backtrajectory Analysis

library(splitr) library(lubridate)

library(devtools)

install_github("rich-iannone/splitr")

Prevent timing out

getOption('timeout') options(timeout=10000)

Build the model

trajectory_model <-hysplit_trajectory( lat = 9.99645,
lon = 76.3427, height = 30, duration = 48, met_type = "reanalysis", direction = "backward", days = seq( lubridate::ymd("2023-01-01"), lubridate::ymd("2023-01-10"), by = "1 day" ), daily_hours = c(8,16,24) )

Remove the missing values (NAs)

trajectory_model_complete <-trajectory_model[complete.cases(trajectory_model), ]

Plot the trajectories

trajectory_plot(trajectory_model_complete)

SGMStalin commented 6 months ago

The same problem.

I solved this by commenting on the following lines in the trajectory_plot() function

(...)

  dates <- traj_df %>% dplyr::pull(traj_dt_i) %>% unique()

  traj_plot <- leaflet::leaflet() %>% 
    leaflet::addProviderTiles(provider = "OpenStreetMap", 
                              group = "OpenStreetMap") %>% 
    leaflet::addProviderTiles(provider = "CartoDB.DarkMatter",
                              group = "CartoDB Dark Matter") %>%
    leaflet::addProviderTiles(provider = "CartoDB.Positron",
                              group = "CartoDB Positron") %>%
    leaflet::addProviderTiles(provider = "Esri.WorldTerrain",
                              group = "ESRI World Terrain") %>%
    # leaflet::addProviderTiles(provider = "Stamen.Toner",
    #                           group = "Stamen Toner") %>%
    leaflet::fitBounds(lng1 = min(traj_df[["lon"]]), 
                       lat1 = min(traj_df[["lat"]]), 
                       lng2 = max(traj_df[["lon"]]), 
                       lat2 = max(traj_df[["lat"]])) %>% 
    leaflet::addLayersControl(baseGroups = c("CartoDB Positron", 
                                             "CartoDB Dark Matter", "Stamen Toner", "ESRI World Terrain"), 
                              overlayGroups = c("trajectory_points", "trajectory_paths"), 
                              position = "topright")

(...)

For that, you need to extract the function, put the cursor on the function (i.e. trajectory_plot(trajectory_model_complete)) and press F2, or CTRL + left click. Copy the function and paste into a new script, assign a name to the function, comment the two lines mentioned above, run the function and prove it.

Ettore91DE commented 4 months ago

Hello everyone,

I am having the exact same issue here. Unfortunately, I am not sure I understood the solution that has been offered..

The same problem.

I solved this by commenting on the following lines in the trajectory_plot() function

(...)

  dates <- traj_df %>% dplyr::pull(traj_dt_i) %>% unique()

  traj_plot <- leaflet::leaflet() %>% 
    leaflet::addProviderTiles(provider = "OpenStreetMap", 
                              group = "OpenStreetMap") %>% 
    leaflet::addProviderTiles(provider = "CartoDB.DarkMatter",
                              group = "CartoDB Dark Matter") %>%
    leaflet::addProviderTiles(provider = "CartoDB.Positron",
                              group = "CartoDB Positron") %>%
    leaflet::addProviderTiles(provider = "Esri.WorldTerrain",
                              group = "ESRI World Terrain") %>%
    # leaflet::addProviderTiles(provider = "Stamen.Toner",
    #                           group = "Stamen Toner") %>%
    leaflet::fitBounds(lng1 = min(traj_df[["lon"]]), 
                       lat1 = min(traj_df[["lat"]]), 
                       lng2 = max(traj_df[["lon"]]), 
                       lat2 = max(traj_df[["lat"]])) %>% 
    leaflet::addLayersControl(baseGroups = c("CartoDB Positron", 
                                             "CartoDB Dark Matter", "Stamen Toner", "ESRI World Terrain"), 
                              overlayGroups = c("trajectory_points", "trajectory_paths"), 
                              position = "topright")

(...)

For that, you need to extract the function, put the cursor on the function (i.e. trajectory_plot(trajectory_model_complete)) and press F2, or CTRL + left click. Copy the function and paste into a new script, assign a name to the function, comment the two lines mentioned above, run the function and prove it.

SGMStalin commented 4 months ago

Can you provide us with more information about what you've done?

Ettore91DE commented 4 months ago

Sure, so I first crated a trajectory model using the hysplit_trajectory function:

`trajectory_model_march2019 <- hysplit_trajectory( lat=51.36588, lon=12.31009, duration=168, height = c(3,15,28,50), met_type="gdas1", days=seq( lubridate::ymd("2019-03-27"), lubridate::ymd("2019-04-02"), by="1 day"), direction="backward", daily_hours=12)

trajectory_model_complete_march2019 <- trajectory_model_march2019[complete.cases(trajectory_model_march2019),] `

then I created a new trajectory plot function: `trajectory_plot_new <- function(x, show_hourly = TRUE, color_scheme = "cycle_hues") {

if (inherits(x, "trajectory_model")) { if (!is.null(x$traj_df)) { traj_df <- x$traj_df } else { stop("There is no data available for plotting.") } }

if (inherits(x, "data.frame")) { if (all(c("run", "receptor", "hour_along", "traj_dt", "lat", "lon", "height", "traj_dt_i") %in% colnames(x))) { traj_df <- x } else { stop("This tibble does not contain plottable trajectory data.") } }

dt_runs <- traj_df$traj_dt_i %>% unique() %>% length()

if (color_scheme == "cycle_hues") { colors <- scales::hue_pal(c = 90, l = 70)(dt_runs) } else if (color_scheme == "increasingly_gray") { colors <- scales::grey_pal(0.7, 0.1)(dt_runs) }

Correct longitude values near prime meridian

traj_df$lon[which(traj_df$lon > 0)] <- traj_df$lon[which(traj_df$lon > 0)] - (180*2)

receptors <- traj_df %>% dplyr::pull(receptor) %>% unique()

dates <- traj_df %>% dplyr::pull(traj_dt_i) %>% unique()

traj_plot <- leaflet::leaflet() %>% leaflet::addProviderTiles( provider = "OpenStreetMap", group = "OpenStreetMap" ) %>% leaflet::addProviderTiles( provider = "CartoDB.DarkMatter", group = "CartoDB Dark Matter" ) %>% leaflet::addProviderTiles( provider = "CartoDB.Positron", group = "CartoDB Positron" ) %>% leaflet::addProviderTiles( provider = "Esri.WorldTerrain", group = "ESRI World Terrain" ) %>%

leaflet::fitBounds(
  lng1 = min(traj_df[["lon"]]),
  lat1 = min(traj_df[["lat"]]),
  lng2 = max(traj_df[["lon"]]),
  lat2 = max(traj_df[["lat"]])
) %>%
leaflet::addLayersControl(
  baseGroups = c(
    "CartoDB Positron", "CartoDB Dark Matter",
   "ESRI World Terrain"
  ),
  overlayGroups = c("trajectory_points", "trajectory_paths"),
  position = "topright"
)

Get different trajectories by receptor and by date

for (i in seq_along(receptors)) {

receptor_i <- receptors[i]

for (j in seq_along(dates)) {

  date_i <- dates[j]

  wind_traj_ij <-
    traj_df %>%
    dplyr::filter(
      receptor == receptor_i,
      traj_dt_i == date_i
    )

  popup <- 
    paste0(
      "<strong>trajectory</strong> ", wind_traj_ij[["traj_dt_i"]],
      "<br><strong>at time</strong> ", wind_traj_ij[["traj_dt"]],
      " (", wind_traj_ij[["hour_along"]],
      " h)<br><strong>height</strong> ", wind_traj_ij[["height"]],
      " <font size=\"1\">m AGL</font> / ",
      "<strong>P</strong> ", wind_traj_ij[["pressure"]],
      " <font size=\"1\">hPa</font>"
    )

  traj_plot <-
    traj_plot %>%
    leaflet::addPolylines(
      lng = wind_traj_ij[["lon"]],
      lat = wind_traj_ij[["lat"]],
      group = "trajectory_paths",
      weight = 2,
      smoothFactor = 1,
      color = colors[j]
    ) %>%
    leaflet::addCircles(
      lng = wind_traj_ij[["lon"]],
      lat = wind_traj_ij[["lat"]],
      group = "trajectory_points",
      radius = 250,
      fill = TRUE,
      color = colors[j],
      fillColor = colors[j], 
      popup = popup
    )
}

}

traj_plot } `

and plotted it: trajectory_plot_new(trajectory_model_complete_march2019)

This way I manage to get the trajectories but these are now split:

is there something in my script that is causing the trajectories to split in the plot?

Thank you very much for your help!

madpk commented 4 months ago

Hi The bug has been identified as Stamen map tiles being moved to Stadia hosting. I have provided an updated R script to circumvent this issue. Check out this video: https://www.youtube.com/watch?v=Y32iJT4bGP4&t=737s

On Mon, May 27, 2024 at 7:41 PM Ettore91DE @.***> wrote:

Sure, so I first crated a trajectory model using the hysplit_trajectory function:

`trajectory_model_march2019 <- hysplit_trajectory( lat=51.36588, lon=12.31009, duration=168, height = c(3,15,28,50), met_type="gdas1", days=seq( lubridate::ymd("2019-03-27"), lubridate::ymd("2019-04-02"), by="1 day"), direction="backward", daily_hours=12)

trajectory_model_complete_march2019 <- trajectory_model_march2019[complete.cases(trajectory_model_march2019),] `

then I created a new trajectory plot function: `trajectory_plot_new <- function(x, show_hourly = TRUE, color_scheme = "cycle_hues") {

if (inherits(x, "trajectory_model")) { if (!is.null(x$traj_df)) { traj_df <- x$traj_df } else { stop("There is no data available for plotting.") } }

if (inherits(x, "data.frame")) { if (all(c("run", "receptor", "hour_along", "traj_dt", "lat", "lon", "height", "traj_dt_i") %in% colnames(x))) { traj_df <- x } else { stop("This tibble does not contain plottable trajectory data.") } }

dt_runs <- traj_df$traj_dt_i %>% unique() %>% length()

if (color_scheme == "cycle_hues") { colors <- scales::hue_pal(c = 90, l = 70)(dt_runs) } else if (color_scheme == "increasingly_gray") { colors <- scales::grey_pal(0.7, 0.1)(dt_runs) } Correct longitude values near prime meridian

traj_df$lon[which(traj_df$lon > 0)] <- traj_df$lon[which(traj_df$lon > 0)] - (180*2)

receptors <- traj_df %>% dplyr::pull(receptor) %>% unique()

dates <- traj_df %>% dplyr::pull(traj_dt_i) %>% unique()

traj_plot <- leaflet::leaflet() %>% leaflet::addProviderTiles( provider = "OpenStreetMap", group = "OpenStreetMap" ) %>% leaflet::addProviderTiles( provider = "CartoDB.DarkMatter", group = "CartoDB Dark Matter" ) %>% leaflet::addProviderTiles( provider = "CartoDB.Positron", group = "CartoDB Positron" ) %>% leaflet::addProviderTiles( provider = "Esri.WorldTerrain", group = "ESRI World Terrain" ) %>%

leaflet::fitBounds( lng1 = min(traj_df[["lon"]]), lat1 = min(traj_df[["lat"]]), lng2 = max(traj_df[["lon"]]), lat2 = max(traj_df[["lat"]]) ) %>% leaflet::addLayersControl( baseGroups = c( "CartoDB Positron", "CartoDB Dark Matter", "ESRI World Terrain" ), overlayGroups = c("trajectory_points", "trajectory_paths"), position = "topright" )

Get different trajectories by receptor and by date

for (i in seq_along(receptors)) {

receptor_i <- receptors[i]

for (j in seq_along(dates)) {

date_i <- dates[j]

wind_traj_ij <- traj_df %>% dplyr::filter( receptor == receptor_i, traj_dt_i == date_i )

popup <- paste0( "trajectory ", wind_traj_ij[["traj_dt_i"]], "
at time ", wind_traj_ij[["traj_dt"]], " (", wind_traj_ij[["hour_along"]], " h)
height ", wind_traj_ij[["height"]], " <font size=\"1\">m AGL / ", "P ", wind_traj_ij[["pressure"]], " <font size=\"1\">hPa" )

traj_plot <- traj_plot %>% leaflet::addPolylines( lng = wind_traj_ij[["lon"]], lat = wind_traj_ij[["lat"]], group = "trajectory_paths", weight = 2, smoothFactor = 1, color = colors[j] ) %>% leaflet::addCircles( lng = wind_traj_ij[["lon"]], lat = wind_traj_ij[["lat"]], group = "trajectory_points", radius = 250, fill = TRUE, color = colors[j], fillColor = colors[j], popup = popup ) }

}

traj_plot } `

and plotted it: trajectory_plot_new(trajectory_model_complete_march2019)

This way I manage to get the trajectories but these are now split: Trajectories.PNG (view on web) https://github.com/rich-iannone/splitr/assets/159437124/455529f2-75f6-45b2-9c86-b8501e5a2049

is there something in my script that is causing the trajectories to split in the plot?

Thank you very much for your help!

— Reply to this email directly, view it on GitHub https://github.com/rich-iannone/splitr/issues/83#issuecomment-2133568150, or unsubscribe https://github.com/notifications/unsubscribe-auth/AVM23Y6TBOZK77V7GAM2HT3ZEM5KRAVCNFSM6AAAAABAEVF23OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMZTGU3DQMJVGA . You are receiving this because you authored the thread.Message ID: @.***>

rich-iannone / splitr

The map appeared but the trajectories not. #83

HYSPLIT Backtrajectory Analysis

library(devtools)

install_github("rich-iannone/splitr")

Prevent timing out

Build the model

Remove the missing values (NAs)

Plot the trajectories

Correct longitude values near prime meridian

Get different trajectories by receptor and by date