R/maps.R

In david6marsh/himach: Find Routes for Supersonic Aircraft

# maps
#
# Functions for plotting nice maps of routes

utils::globalVariables(c("crs_longlat"))
utils::globalVariables(c("mach_kph"))

#' Version of \code{st_transform} with view window to avoid dateline
#'
#' \code{st_window} does a \code{st_transform} but first cuts the data to an
#' appropriate view window and so avoids problems with objects wrapping around
#' the back of the globe
#'
#'
#' \code{\link[sf:st_transform]{st_wrap_dateline}} _should_ handle the break in
#' a map projections but uses `GDAL` for this. Given persistent issues in
#' installing GDAL, \code{st_window} achieves the same using \code{s2} instead.
#'
#' It works for any 'simple' projection, in the sense of one that has a dateline
#' that is a single line of longitude: ie the proj4string contains either
#' "longitude_of_center", so the dateline is that +180; or not, in which case it
#' assumes the "longitude_of_center" is 0.
#'
#'
#' @param m A map dataframe, ie of class \code{sf} and \code{data.frame}, or an
#'   \code{sfc_MULTIPOLYGON}
#' @param crs Destination coordinate reference system, as in \code{st_tranform}
#' @param longit_margin Amount trimmed off the 'far side' of the projection in
#'   degrees.
#'
#' @return \code{sf} dataframe, same as the parameter \code{m}
#'
#' @import sf
#' @import dplyr
#' @import tidyr
#'
#' @examples
#' world <- sf::st_as_sf(rnaturalearthdata::coastline110)
#' w_pacific <- st_window(world, crs_Pacific)
#' ggplot2::ggplot(w_pacific) + ggplot2::geom_sf()
#'
#' # bad - not run - dateline problem example
#' # ggplot2::ggplot(st_transform(world, crs_Pacific)) +
#' #   ggplot2::geom_sf()
#'
#' @export
st_window <- function(m, crs = himach::crs_Atlantic, longit_margin = 0.1){
  # 'dateline' of the crs
  longit <- mod_long(long_cent(crs) + 180)
  long1 <- mod_long(longit + longit_margin)
  long2 <- mod_long(longit - longit_margin)
  #make window: oriented = TRUE because otherwise so big, s2 assumes inverse
  window <- s2::s2_make_polygon(c(rep(long1, 181),
                                  rep(long2, 181)),
                                c(seq.int(90, -90, -1),  seq.int(-90, 90, 1)),
                                oriented = TRUE)
  m %>%
    sf::st_as_s2() %>%
    s2::s2_intersection(window,
                        s2::s2_options(model = "open")) %>%
    sf::st_as_sfc() %>%
    sf::st_transform(crs)
}

#wrapper around ggplot/geom_sf for quick map
plot_map <- function(msf,
                     c_land = "green1",
                     c_border = "grey70",
                     w_border = 0.1){
  ggplot2::ggplot(msf) +
    ggplot2::geom_sf(linewidth = w_border, colour = c_border, fill = c_land) +
    # coord_sf() + #not needed for already-projected feature sets
    ggplot2::theme_minimal()
}


# achievable range from a point
# not used for route finding
#' @import sf
make_range_envelope <- function(ac, ap, ap_locs = make_airports(),
                                envelope_points=70){
  #range envelope shows how far from an airport you can go  with a given range
  #create no-wind range ellipse

  ap_loc <- ap_locs %>%
    filter(.data$APICAO == ap)
  # centre of range
  geo_c <- c(ap_loc$long, ap_loc$lat)

  # use CRS centred on centre of route envelopes
  cen_prj <- sf::st_crs(paste0("+proj=laea +lat_0=", round(geo_c[2],1),
                               " +lon_0=", round(geo_c[1],1),
                               " +x_0=4321000 +y_0=3210000 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
  dist <- ac$range_km * 1000
  theta <- seq(0, 360, length.out = envelope_points)

  geod <- geosphere::geodesic(geo_c, theta, dist)

  # convert to simple feature
  pg <- sf::st_multipoint(geod[ ,1:2]) %>%
    sf::st_cast('LINESTRING') %>%
    sf::st_cast('POLYGON') %>%
    sf::st_sfc(crs = crs_longlat) %>%
    sf::st_transform(cen_prj) %>%
    # occasionally fails as self-intersection when later st_intersection
    # so this should solve that
    sf::st_make_valid()
}

# simplify a map
thin <- function(m, tol_km = 4){
  m  %>%
    st_as_s2() %>%
    s2::s2_simplify(tolerance = tol_km * 1000) %>%
    st_as_sfc()
}

#' Map a set of routes
#'
#' \code{map_routes} plots routes, with many options
#'
#' This function plots the routes, with options for additional layers. Multiple
#' routes are expected, and they can be coloured by time advantage, by speed
#' along each segment, or by aircraft type.
#'
#' The option \code{show_route} "time" requires 'advantage_h' to have been added
#' to the routes set, from the route summary. If it hasn't then this is done in
#' a local version, then discarded. Running \code{summarise_routes} to do this
#' requires an airport dataset; if \code{is.na(ap_loc)} then this is not
#' available, so a default set is used. You can turn on \code{warn} to see if
#' this is happening, but by default it is silent.
#'
#' @param thin_map The minimum is a \code{MULTIPOLYGON} map, 'thin' in that it
#'   is without buffer, so a normal coastline map.
#' @param routes as generated by \code{\link{find_route}}
#' @param crs Coordinate reference system, default \code{crs_Atlantic}.
#' @param show_route  one of "speed", "aircraft", "time", "circuity", "accel",
#'   "traffic" to indicate what goes in the legend.
#' @param fat_map optional coast + buffer map, default NA.
#' @param forecast,fc_var,fc_text Forecast set and two strings. See details,
#'   default NA.
#' @param avoid_map optional map of no-fly zones, default NA.
#' @param ap_loc Show used origin and destination airports if this is a set of
#'   airports from \code{\link{make_airports}}, or not if NA (default). This
#'   dataset can be all airports, and is filtered to those used by
#'   \code{routes}.
#' @param ap_col,ap_size Colour and size of used airport markers (dark blue,
#'   0.4)
#' @param refuel_airports Show the used refuel airports using these locations,
#'   or nothing if NA. (Defaults to same as \code{ap_loc}.)
#' @param rap_col,rap_size Colour and size of refuel airport markers (red, 0.4)
#' @param crow,crow_col,crow_size If TRUE, show the 'crow-flies' direct great
#'   circle, in colour \code{crow_col} and thickness \code{crow_size}. Default
#'   FALSE, "grey70", 0.2
#' @param route_envelope show the route envelope (default FALSE).
#' @param e_col,e_alpha,e_size colour, alpha and width for the range envelope.
#'   Default "grey70", 0.4, 0.6
#' @param bound,bound_margin_km If bound=TRUE (default) crop to bounding box of
#'   the \code{routes}, with additional \code{bound_margin_km} in km (default
#'   200)
#' @param simplify_km Simplify the two maps to this scale before plotting
#'   (default 10).
#' @param land_f,buffer_f,avoid_f fill colours for thin, fat and no-fly maps,
#'   default grey 90, 70 and 80, respectively
#' @param land_c,land_s boundary colour and size for land areas (countries),
#'   default grey 85 and 0.2, respectively (use NA to turn off)
#' @param avoid_c,avoid_s boundary colour and size for avoid areas, default grey
#'   95 and 0.3, respectively
#' @param l_alpha,l_size line (route) settings for alpha (transparency) and
#'   width, defaults 0.6 and 0.4.
#' @param scale_direction Passed to scale_colour_viridis, either -1 (default) or
#'   or 1.
#' @param title,subtitle Passed to ggplot.
#' @param warn if TRUE show some warnings (when defaults loaded) (default FALSE)
#' @param ... further parameters passed to \code{scale_colour_viridis_b} (or _c,
#'   _d), such as \code{breaks = }.
#'
#' @return A \code{ggplot}.
#'
#' @details
#'
#' For \code{show_route = } "speed", "aircraft", "time", "circuity" or "accel",
#' the information is already available in the \code{routes} dataset. For
#' \code{show_route = "traffic"} you need to provide a forecast dataset that
#' contains at least the \code{fullRouteID and acID} fields which are normal in
#' the \code{routes} dataset, and a field giving the volume of the forecast
#' \code{fc_var}. This could be flights, seats, or something else: use
#' \code{fc_text} for the legend title to show the units of \code{fc_var}.
#' Combinations of \code{fullRouteID and acID} must be unique, which probably
#' means you must filter by forecast year and forecast scenario before passing
#' to \code{map_routes}.
#'
#' The time to compute the map may not be very different with \code{simplify_km}
#' varying between 2km and 20km, but the time to plot on the screen, or
#' \code{ggsave} to a file, is longer than the compute time. It is this latter
#' time that's reduced by simplifying the maps. For single, or short routes, you
#' can probably see the difference between 2km and 10km, so it's your choice to
#' prefer speed or beauty.
#'
#' @import sf
#' @import dplyr
#' @import tidyr
#' @import ggplot2
#'
#' @examples
#' #see introductory vignette
#'
#' @export
map_routes <- function(
    thin_map, routes=NA, crs=himach::crs_Atlantic,
    show_route = c("speed", "aircraft","time",
                   "circuity", "acceleration", "traffic"),
    fat_map=NA, avoid_map=NA,
    ap_loc=NA, ap_col="darkblue", ap_size=0.4,
    forecast=NA, fc_var = NA_character_, fc_text = NA_character_,
    crow=FALSE, crow_col="grey70", crow_size=0.2,
    route_envelope=FALSE,
    bound=TRUE, bound_margin_km=200,
    simplify_km = 8,
    land_f="grey90", buffer_f="grey60",
    land_c="grey85", land_s=0.2,
    avoid_f="grey80", avoid_c="grey95", avoid_s=0.3,
    l_alpha=0.8, l_size=0.5,
    e_alpha=0.4, e_size=0.6, e_col="grey70",
    refuel_airports=ap_loc, rap_col="red", rap_size=0.4,
    scale_direction = -1,
    title = "", subtitle = "",
    warn = FALSE,
    ...
){
  # use tidyverse approach https://design.tidyverse.org/def-enum.html
  show_route <- match.arg(show_route)

  # remove the non-routes (have time = NA)
  # these are where refuelling was needed
  if (is.data.frame(routes)) routes <- routes %>% filter(!is.na(time_h))

  # older sets of routes may not be sf
  if (is.data.frame(routes) && !"sf" %in% class(routes)){
    if (warn) message("routes must be class sf. Doing this for you.")
    routes <- st_set_geometry(routes, "gc")
  }

  #thin map is the one without buffer
  # thin_map <- st_window(thin_map, crs) #force CRS and cut to view window

  #layer 1 (one or two base maps)
  if (is.na(fat_map)) {
    m <- plot_map(thin_map %>% thin(simplify_km) %>% st_window(crs),
                  c_border=land_c, c_land=land_f, w_border = land_s)
  } else {
    m <- plot_map(fat_map %>% thin(simplify_km) %>% st_window(crs),
                  c_border=NA, c_land=buffer_f) +
      geom_sf(data=thin_map %>% thin(simplify_km) %>% st_window(crs), fill=land_f,
              colour=land_c, linewidth = land_s)
  }

  #layer 2 (no fly-zone)
  if (is.list(avoid_map)){
    m <- m +
      geom_sf(data = st_window(avoid_map, crs),
              colour=avoid_c, fill=avoid_f, linewidth = avoid_s)
  }

  # 3: prelim: check if summarise_routes has been run
  # by seeing if 'advantage_h' has been calculated
  if (is.data.frame(routes)) {
    #remove empty routes
    routes <- routes %>%
      filter(!is.na(time_h))

    if ( !("advantage_h" %in% names(routes))) {
      if (warn) message("Adding advantage_h temporarily to the routes set.")
      airports <- ap_loc
      if (!is.data.frame(airports)){
        airports <- make_airports(warn = FALSE)
        if (warn) message("Using default airport set for temporary route summary.")
      }
      rtes <- summarise_routes(routes, airports) %>%
        st_set_geometry(NULL) %>% # for left_join, must _not_ be sf
        select("fullRouteID", "advantage_h", "circuity")

      routes <- routes %>%
        left_join(rtes, by = "fullRouteID") %>%
        arrange(.data$advantage_h)

      routes$gc <- st_window(routes$gc, crs)
      routes <- st_set_geometry(routes, "gc") # default geometry for plots
    }

    #layer 3 (main lines)
    if (show_route=="time"){
      m <- m +  labs(colour = "Time Advantage") +
        geom_sf(data = routes,
                aes(colour = .data$advantage_h),
                fill="white",
                linewidth=l_size, lineend="round", alpha=l_alpha) +
        scale_colour_viridis_c(direction = scale_direction, ...)
    }
    if (show_route=="circuity"){
      m <- m +  labs(colour = "Circuity\n(best=0)") +
        geom_sf(data = routes,
                aes(colour = .data$circuity), fill = "white",
                linewidth=l_size, lineend="round", alpha=l_alpha) +
        scale_colour_viridis_c(direction = scale_direction,
                               labels = scales::percent, ...)
    }
    if (show_route == "speed"){
      m <- m +  labs(colour = "Average speed on segment (kph)") +
        geom_sf(data = routes,
                aes(colour = .data$speed_kph),
                linewidth=l_size, lineend="round", alpha=l_alpha) +
        scale_colour_viridis_c(direction = scale_direction, ...)
    }
    if (show_route == "aircraft"){
      m <- m +  labs(colour = "Aircraft") +
        geom_sf(data = routes,
                aes(colour = .data$acID),
                linewidth=l_size, lineend="round", alpha=l_alpha,
                show.legend = "line")+
        scale_colour_viridis_d(direction = scale_direction, ...)
    }
    if (show_route == "accel"){
      m <- m +  labs(colour = "Number of accelerations\nto supersonic") +
        geom_sf(data = routes,
                aes(colour = factor(.data$n_accel)),
                linewidth=l_size, lineend="round", alpha=l_alpha,
                show.legend = "line") +
        scale_colour_viridis_d(direction = scale_direction, ...)
    }
    if (show_route == "traffic"){
      # check forecast has right format
      stopifnot (is.data.frame(forecast), !is.na(fc_var), !is.na(fc_text)) # must have provided forecast
      stopifnot("fullRouteID" %in% names(forecast), "acID" %in% names(forecast), fc_var %in% names(forecast))
      if (nrow(forecast %>% group_by(.data$fullRouteID, .data$acID) %>%
               summarise(n = n()) %>% filter(n > 1)) > 0) {
        stop("Rows in forecast not unique per route & aircraft. Have you filtered on year and scenario?")
      }
      fc <- forecast %>%
        ungroup() %>%
        select("fullRouteID", "acID", {{ fc_var }})
      # check if routes missing for the forecast
      route_miss <- fc %>%
        anti_join(routes %>% st_drop_geometry() %>% select("fullRouteID", "acID"),
                  by = c("fullRouteID", "acID"))
      if (nrow(route_miss) > 0) {
        warning(nrow(route_miss), " forecast rows do not have routes available.", call. = FALSE)
      }
      m <- m +  labs(colour = fc_text) +
        geom_sf(data = routes %>%
                  inner_join(fc, by = c("fullRouteID", "acID")) %>%
                  arrange(.data[[fc_var]]), # plot in increasing order
                aes(colour = .data[[fc_var]]),
                linewidth=l_size, lineend="round", alpha=l_alpha) +
        scale_colour_viridis_b(direction = scale_direction, ...)
    }
  }


  #layer 4: crow-flies
  if (crow){
    m <- m +
      geom_sf(data = st_window(routes$crow, crs),
              colour=crow_col, linewidth = crow_size)
  }

  #layer 5: range envelope
  if (route_envelope){
    m <- m +
      geom_sf(data = st_window(st_cast(routes$envelope, 'MULTILINESTRING'),
                               crs),
              fill = NA,
              colour = e_col, alpha = e_alpha, linewidth = e_size)
  }

  #layer 6: airports
  if (is.data.frame(ap_loc)){
    used_APs <- sort(unique(unlist(lapply(unique(routes$routeID),
                                          function(x)strsplit(x, "<>")))))
    APs <- ap_loc %>%
      filter(.data$APICAO %in% used_APs)
    APs$ap_locs <- prj(APs$ap_locs, crs = crs)
    APs <- st_set_geometry(APs, "ap_locs")
    m <- m +
      geom_sf(data = APs,
              colour=ap_col, size=ap_size)
  }

  #layer 7: refuel airports
  if (is.data.frame(refuel_airports)){
    used_refuel_APs <- sort(unique(routes$refuel_ap))
    rAPs <- refuel_airports %>%
      filter(.data$APICAO %in% used_refuel_APs)
    rAPs$ap_locs <- prj(rAPs$ap_locs, crs = crs)
    rAPs <- st_set_geometry(rAPs, "ap_locs")
    m <- m +
      geom_sf(data = rAPs,
              colour=rap_col, size=rap_size)
  }

  #apply bounds?
  if (is.data.frame(routes) && bound){
    #crop based on the bounding box of all of the routes
    bbox <- st_bbox(st_transform(routes$gc, crs=crs)) +
      c(-1,-1,1,1) * bound_margin_km * 10^3 #zoomed box + margin-km
    m <- m +
      xlim(bbox$xmin, bbox$xmax) + ylim(bbox$ymin, bbox$ymax)
  }

  m <- m +
    labs(title = title, subtitle = subtitle) +
    theme(legend.position = "bottom")

  return(m)
}

#find the longitude of centre
long_cent <- function(m){
  #full text version of crs
  tx <- st_as_text(st_crs(m))
  #search as vector
  txv <- t(stringr::str_split(tx,",", simplify = TRUE))
  txv <- stringr::str_remove_all(txv, "\\]") #get rid of these
  # centre is made up of primemeridiand and long of center if either exists
  w <- which(stringr::str_detect(txv, "PRIMEM"))
  if (length(w) == 0){
    # if not mentioned, default to 0
    pm <- 0
  } else {
    pm <- as.numeric(txv[w + 1])
  }

  w <- which(stringr::str_detect(txv, "longitude_of_center"))
  if (length(w) == 0){
    # if not mentioned, default to 0
    c_long <- 0
  } else {
    c_long <- as.numeric(txv[w + 1])
  }
  pm + c_long
}


#' Profile a set of routes
#'
#' @param routes as generated by \code{\link{find_route}}
#' @param yvar horizontal axis is hours or longitude
#' @param ap_loc Airports and coordinates, by (silent) default from \code{\link{make_airports}}
#' @param n_max maximum number of routes to plot (default 2)
#'
#' @return A list of named list pairs of plots, which can be displayed using eg \code{result[1]}.
#'
#' @import ggplot2
#' @import dplyr
#'
#' @examples
#' # not run ---
#' # plot_list <- profile_routes(routes, n_max = 3)
#' # plot_list # to display them all
#'
#' @export
profile_routes <- function(routes, yvar = c("hours", "longitude"),
                           ap_loc = make_airports(warn = FALSE), n_max = 2){
  yvar <- match.arg(yvar) # check the yvar argument
  xtitle <- ifelse(yvar == "hours",
                   "Elapsed time (hours)",
                   "Longitude (degrees)")
  routes <- routes %>%
    filter(!is.na(.data$time_h))
  # check reasonable number of plots
  route_id <- unique(routes$routeID)
  if (length(route_id) == 0) stop("No flyable routes to profile.")
  if (length(route_id) > n_max) {
    warning("Too many routes. Just profiling n_max=", n_max, ".")
    route_id <- route_id[1:n_max]
  }

  rr <- routes |>
    filter(!is.na(.data$time_h) & .data$routeID %in% route_id) |>
    group_by(.data$timestamp) |>
    # get airport labels - this seems ridiculous
    filter(!is.na(.data$time_h)) |>
    left_join(ap_loc |> select(-"ap_locs"), by = c("from_long"="long", "from_lat"="lat")) |>
    rename(ap1_name = "APICAO") |>
    left_join(ap_loc |> select(-"ap_locs"), by = c("to_long"="long", "to_lat"="lat")) |>
    rename(ap2_name = "APICAO") |>
    mutate(lab = coalesce(.data$ap1_name, .data$ap2_name)) |>
    select(-"ap1_name", -"ap2_name") |>
    mutate(lab = case_when(
      row_number() == 1 ~ .data$lab,
      lag(.data$lab) == .data$lab ~ "",
      is.na(.data$lab) ~ "",
      TRUE ~ .data$lab)) |>
    # calculate plot points
    mutate(speed_M = .data$speed_kph/mach_kph,
           elapsed_h = cumsum(.data$time_h),
           yvar = yvar,
           xvalue = if_else(yvar=="hours",
                            coalesce(lag(.data$elapsed_h),0),
                            .data$from_long),
           xvalue2 = if_else(yvar=="hours",
                             .data$elapsed_h,
                             .data$to_long),
           cum_dist_km = cumsum(.data$gcdist_km),
           prev_dist = coalesce(lag(.data$cum_dist_km),0)) |>
    # label position varies
    mutate(lab_h = case_when(
      row_number() == 1 ~ .data$xvalue,
      TRUE ~ .data$xvalue2),
      lab_d = case_when(
        row_number() == 1 ~ .data$xvalue,
        TRUE ~ .data$cum_dist_km))

  lapply(route_id, function(rte){
    rr1 <- rr |> filter(.data$routeID == rte)
    # make label
    ap_pair <- stringr::str_c(rr1$lab[1], "<>", rr1$lab[nrow(rr1)])
    #distance v time
    (dvt <- ggplot(rr1, aes(x=.data$xvalue, xend=.data$xvalue2,
                            y=.data$prev_dist, yend=.data$cum_dist_km, colour=.data$acID)) +
        geom_segment(linewidth=1) +
        labs(x=xtitle, y="Flown distance (km)",
             colour="Aircraft") +
        theme_minimal() + theme(legend.position="top"))
    #speed v time
    (svt <- ggplot(rr1, aes(x=.data$xvalue, xend=.data$xvalue2,
                            y=.data$speed_M, yend=.data$speed_M, colour=.data$acID)) +
        geom_hline(yintercept = 1.0, colour="grey70") + #supersonic boundary
        geom_segment(linewidth=1) +
        geom_text(aes(x=.data$lab_h, y=0, label=.data$lab)) +
        labs(y="Ave. speed (Mach)", x = NULL,
             colour="Aircraft",
             subtitle = stringr::str_c(ap_pair, " Performance v. Time"),
             caption=NULL) +
        scale_y_continuous(limits=c(0,NA)) + #include 0
        theme_minimal() + guides(colour="none"))
    #speed v distance
    (svd <- ggplot(rr1, aes(y=.data$speed_M, yend=.data$speed_M,
                            x=.data$prev_dist, xend=.data$cum_dist_km, colour=.data$acID)) +
        geom_hline(yintercept = 1.0, colour="grey70") + #supersonic boundary
        geom_segment(linewidth=1) +
        geom_text(aes(x=.data$lab_d, y=0, label=.data$lab)) +
        labs(y="Ave. speed (Mach)", x = NULL,
             colour="Aircraft",
             subtitle = stringr::str_c(ap_pair, " Performance v. Distance Flown"),
             caption=NULL) +
        scale_y_continuous(limits=c(0,NA)) + #include 0
        theme_minimal() + guides(colour="none"))

    list(time = cowplot::plot_grid(svt, dvt, align=c("v"), ncol=1),
         distance = cowplot::plot_grid(svd, dvt+coord_flip(), align=c("v"), ncol = 1))
  })
}