4

enter image description hereI have been meaning to create a new geom for a data set that has been tidied in the following form:

Katrina
# A tibble: 3 x 9
  storm_id     date                latitude longitude wind_speed    ne    se    sw    nw
  <chr>        <dttm>                 <dbl>     <dbl> <fct>      <dbl> <dbl> <dbl> <dbl>
1 KATRINA-2005 2005-08-29 12:00:00     29.5     -89.6 34           200   200   150   100
2 KATRINA-2005 2005-08-29 12:00:00     29.5     -89.6 50           120   120    75    75
3 KATRINA-2005 2005-08-29 12:00:00     29.5     -89.6 64            90    90    60    60

I first defined the class and then the actual geom function, however, my output plot turns out to be so miniaturized, So I would appreciate it if you could tell me where I possibly go wrong with the scales.

GeomHurricane <- ggplot2::ggproto("GeomHurricane", Geom, 
                         required_aes = c("x", "y",
                                          "r_ne", "r_se", "r_sw", "r_nw"
                                          ),
                         default_aes = aes(fill = 1, colour = 1, 
                                           alpha = 1, scale_radii = 1),
                         draw_key = draw_key_polygon,
                         
                         draw_group = function(data, panel_scales, coord) {
                           
                            coords <- coord$transform(data, panel_scales) %>%
                              mutate(r_ne = r_ne * 1852 * scale_radii, 
                                     r_se = r_se * 1852 * scale_radii, 
                                     r_sw = r_sw * 1852 * scale_radii,
                                     r_nw = r_nw * 1852 * scale_radii
                                                   )
                           
                           # Creating quadrants 
                           for(i in 1:nrow(data)) {
                             
                             # Creating the northeast quadrants
                             data_ne <- data.frame(colour = data[i,]$colour,
                               fill = data[i,]$fill,
                               geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                    b = 1:90,
                                                    d = data[i,]$r_ne),
                               group = data[i,]$group,
                               PANEL = data[i,]$PANEL,
                               alpha = data[i,]$alpha
                             )
                             
                             # Creating the southeast quadrants
                             data_se <- data.frame(colour = data[i,]$colour, 
                               fill = data[i,]$fill,
                               geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                    b = 90:180,
                                                    d = data[i,]$r_se),
                               group = data[i,]$group,
                               PANEL = data[i,]$PANEL,
                               alpha = data[i,]$alpha
                             )
                             
                             # Creating the southwest quadrants
                             data_sw <- data.frame(colour = data[i,]$colour, 
                               fill = data[i,]$fill,
                               geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                    b = 180:270,
                                                    d = data[i,]$r_sw),
                               group = data[i,]$group,
                               PANEL = data[i,]$PANEL,
                               alpha = data[i,]$alpha
                             )
                             
                             # Creating the northwest quadrants
                             data_nw <- data.frame(colour = data[i,]$colour,
                               fill = data[i,]$fill, 
                               geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                    b = 270:360,
                                                    d = data[i,]$r_nw),
                               group = data[i,]$group,
                               PANEL = data[i,]$PANEL,
                               alpha = data[i,]$alpha
                             )
                             
                             data_quadrants <- dplyr::bind_rows(list(
                               data_ne, data_se, data_sw, data_nw
                             )) 
                             
                             data_quadrants <- data_quadrants %>% dplyr::rename(
                               x = lon,
                               y = lat
                             )
                             
                             data_quadrants$colour <- as.character(data_quadrants$colour)
                             data_quadrants$fill <- as.character(data_quadrants$fill)
                             
                           }

                             coords_data <- coord$transform(data_quadrants, panel_scales)
                             
                             grid::polygonGrob(
                               x = coords_data$x,
                               y = coords_data$y,
                               default.units = "native", 
                               gp = grid::gpar(
                                 col = coords_data$colour, 
                                 fill = coords_data$fill,
                                 alpha = coords_data$alpha
                               )
                             )
                          }
)

and the actual geom function definition:

geom_hurricane <- function(mapping = NULL, data = NULL, stat = "identity",
                           position = "identity", na.rm = FALSE,
                           show.legend = NA, inherit.aes = TRUE, ...) {
  ggplot2::layer(
    geom = GeomHurricane, mapping = mapping,
    data = data, stat = stat, position = position,
    show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}

So I went on to plot the following:

ggplot(data = Katrina) + 
  geom_hurricane(aes(x = longitude, y = latitude, 
                     r_ne = ne, r_se = se, r_sw = sw, r_nw = nw,
                     fill = wind_speed, colour = wind_speed)) + 
  scale_colour_manual(name = "Wind speed (kts)",
                      values = c("red", "orange", "yellow")) +
  scale_fill_manual(name = "Wind speed (kts)",
                    values = c("red", "orange", "yellow"))

The data for this purpose can be found here as Atlantic basin data set 1988 - 2018: http://rammb.cira.colostate.edu/research/tropical_cyclones/tc_extended_best_track_dataset/

For your consideration I used the following codes to tidy the data:

ext_tracks_widths <- c(7, 10, 2, 2, 3, 5, 5, 6, 4, 5, 4, 4, 5, 3, 4, 3, 3, 3,
                       4, 3, 3, 3, 4, 3, 3, 3, 2, 6, 1)


ext_tracks_colnames <- c("storm_id", "storm_name", "month", "day",
                         "hour", "year", "latitude", "longitude",
                         "max_wind", "min_pressure", "rad_max_wind",
                         "eye_diameter", "pressure_1", "pressure_2",
                         paste("radius_34", c("ne", "se", "sw", "nw"), sep = "_"),
                         paste("radius_50", c("ne", "se", "sw", "nw"), sep = "_"),
                         paste("radius_64", c("ne", "se", "sw", "nw"), sep = "_"),
                         "storm_type", "distance_to_land", "final")

ext_tracks <- read_fwf("ebtrk_atlc_1988_2015.txt",
                       fwf_widths(ext_tracks_widths, ext_tracks_colnames), 
                       na = "-99")

storm_observation <- ext_tracks %>%
  unite("storm_id", c("storm_name", "year"), sep = "-", 
        na.rm = TRUE, remove = FALSE) %>%
  mutate(longitude = -longitude) %>%
  unite(date, year, month, day, hour) %>%
  mutate(date = ymd_h(date)) %>%
  select(storm_id, date, latitude, longitude, radius_34_ne:radius_64_nw) %>%
  pivot_longer(cols = contains("radius"), names_to = "wind_speed", 
               values_to = "value") %>%
  separate(wind_speed, c(NA, "wind_speed", "direction"), sep = "_") %>%
  pivot_wider(names_from = "direction", values_from = "value") %>%
  mutate(wind_speed = as.factor(wind_speed))


Katrina <- storm_observation %>%
  filter(storm_id == "KATRINA-2005", date == ymd_h("2005-08-29-12"))
Anoushiravan R
  • 21,622
  • 3
  • 18
  • 41
  • Hey this sounds like an interesting question. Would you mind including a picture of the output and describing how it should be changed? – teunbrand Jan 17 '21 at 14:18
  • Yes, Sure. Thank you very much. If there is any clarification I need to make feel free to ask. By the way I am doing this assignment as part of my R self-study learning, So I would like to fully grasp the concept since the documentation about building geom is not very clear and helpful on some aspects of it, I became a bit desperate. – Anoushiravan R Jan 17 '21 at 14:41
  • Yeah I know the documentation is a bit sparse. Here are two places that were useful to me: https://ggplot2.tidyverse.org/articles/extending-ggplot2.html and https://ggplot2-book.org/extensions.html. I'll have a look! – teunbrand Jan 17 '21 at 14:45
  • Thank you very much, that's very kind of you. – Anoushiravan R Jan 17 '21 at 15:23

1 Answers1

3

Alright, there are 2 problems that I spotted. Problem 1 is that in your draw_group() ggproto method, you convert the radii from nautical miles to meters (I think), but you write this to the coords variable. However, you use the data variable to do the geosphere::destPoint calculation.

Here is a version of that method that I think should work:

  draw_group = function(data, panel_scales, coord) {

    scale_radii <- if (is.null(data$scale_radii)) 1 else data$scale_radii
    data <- data %>%
      mutate(r_ne = r_ne * 1852 * scale_radii, 
             r_se = r_se * 1852 * scale_radii, 
             r_sw = r_sw * 1852 * scale_radii,
             r_nw = r_nw * 1852 * scale_radii
      )
    
    # Creating quadrants 
    for(i in 1:nrow(data)) {
      
      # Creating the northeast quadrants
      data_ne <- data.frame(colour = data[i,]$colour,
                            fill = data[i,]$fill,
                            geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                 b = 1:90, # Should this start at 0?
                                                 d = data[i,]$r_ne),
                            group = data[i,]$group,
                            PANEL = data[i,]$PANEL,
                            alpha = data[i,]$alpha
      )
      
      # Creating the southeast quadrants
      data_se <- data.frame(colour = data[i,]$colour, 
                            fill = data[i,]$fill,
                            geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                 b = 90:180,
                                                 d = data[i,]$r_se),
                            group = data[i,]$group,
                            PANEL = data[i,]$PANEL,
                            alpha = data[i,]$alpha
      )
      
      # Creating the southwest quadrants
      data_sw <- data.frame(colour = data[i,]$colour, 
                            fill = data[i,]$fill,
                            geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                 b = 180:270,
                                                 d = data[i,]$r_sw),
                            group = data[i,]$group,
                            PANEL = data[i,]$PANEL,
                            alpha = data[i,]$alpha
      )
      
      # Creating the northwest quadrants
      data_nw <- data.frame(colour = data[i,]$colour,
                            fill = data[i,]$fill, 
                            geosphere::destPoint(p = c(data[i,]$x, data[i,]$y),
                                                 b = 270:360,
                                                 d = data[i,]$r_nw),
                            group = data[i,]$group,
                            PANEL = data[i,]$PANEL,
                            alpha = data[i,]$alpha
      )
      
      data_quadrants <- dplyr::bind_rows(list(
        data_ne, data_se, data_sw, data_nw
      )) 
      
      data_quadrants <- data_quadrants %>% dplyr::rename(
        x = lon,
        y = lat
      )
      
      data_quadrants$colour <- as.character(data_quadrants$colour)
      data_quadrants$fill <- as.character(data_quadrants$fill)
      
    }
    
    coords_data <- coord$transform(data_quadrants, panel_scales)
    
    grid::polygonGrob(
      x = coords_data$x,
      y = coords_data$y,
      default.units = "native", 
      gp = grid::gpar(
        col = coords_data$colour, 
        fill = coords_data$fill,
        alpha = coords_data$alpha
      )
    )
  }

The next problem is that you only define 1 x coordinate with the Katrina example. However, the scales don't know about your radius parameters, so they don't adjust the limits to fit your radii in. You can circumvent this by setting xmin, xmax, ymin and ymax bounding box parameters, so that scale_x_continuous() can learn about your radius. (Same thing for the y scale). I'd solve this by using a setup_data method to your ggproto object.

Here is the setup data method that I used to test with, but I'm no spatial genius so you'd have to check if this makes sense.

  setup_data = function(data, params) {

    maxrad <- max(c(data$r_ne, data$r_se, data$r_sw, data$r_nw))
    maxrad <- maxrad * 1852

    x_range <- unique(range(data$x))
    y_range <- unique(range(data$y))
    xy <- as.matrix(expand.grid(x_range, y_range))

    extend <- lapply(c(0, 90, 180, 270), function(b) {
      geosphere::destPoint(p = xy,
                           b = b,
                           d = maxrad)
    })
    extend <- do.call(rbind, extend)

    transform(
      data,
      xmin = min(extend[, 1]),
      xmax = max(extend[, 1]),
      ymin = min(extend[, 2]),
      ymax = max(extend[, 2])
    )
  }

After implenting these changes, I get a figure like this:

enter image description here

teunbrand
  • 33,645
  • 4
  • 37
  • 63
  • Thank you very much dear @teunbrand, Do you suggest I should define a new stat and put your `setup_data` in it? – Anoushiravan R Jan 17 '21 at 16:41
  • You you can just define the `setup_data` in the ggproto of the geom itself if you want: it is part of every geom. Just print `GeomPoint$setup_data` at the console for example and you'll see that `geom_point()` also has one (but doesn't do anything useful). – teunbrand Jan 17 '21 at 16:43
  • Oh man! I don't know how to thank you! I honestly think you are indeed a spatial genius. I got what I have been looking for however all of these are a bit too much for me as I have been learning R for almost 7 months. I would be grateful If you could kindly tell me what I need to do improve my grasp on building geom subject. – Anoushiravan R Jan 17 '21 at 16:58
  • This is my fourth coursera course and the only one I struggled and didn't know what I should do and come up with my own solution. So I sought help. However, I guess I am expecting too much of myself. I would like to learn day by day and look up to someone like you as inspiration. – Anoushiravan R Jan 17 '21 at 17:00
  • I don't know many people who write their own geoms and ggproto questions are not very common on SO. The tricky thing with writing your own geoms is to debug problems that you might have and I gave a few tips over [here](https://stackoverflow.com/a/62702645/11374827) that you might find convenient. – teunbrand Jan 17 '21 at 17:30
  • Also, the programming style of ggproto is very close to R6 object oriented programming, which you can read about here: https://adv-r.hadley.nz/r6.html – teunbrand Jan 17 '21 at 17:31
  • Thank you very much indeed for your fantastic comments. I also plan on reading Hadley's books ggplot2 and advanced R. I guess They are essentials. – Anoushiravan R Jan 17 '21 at 17:55