建立一个新的 geom_hurricane
Building a new geom_hurricane
我一直想为整理后的数据集创建一个新的geom，形式如下：
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
我首先定义了 class 然后是实际的 geom 函数，然而，我的输出图变得如此小型化，所以如果你能告诉我我可能哪里出错了，我将不胜感激体重秤。
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
)
)
}
)
和实际几何函数定义：
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, ...)
)
}
所以我继续绘制以下内容：
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"))
用于此目的的数据可在此处找到，作为大西洋盆地数据集 1988 - 2018：
http://rammb.cira.colostate.edu/research/tropical_cyclones/tc_extended_best_track_dataset/
为了您的考虑，我使用了以下代码来整理数据：
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"))
好的，我发现了 2 个问题。问题 1 是，在您的 draw_group() ggproto 方法中，您将半径从海里转换为米（我认为），但是您将其写入 coords 变量。但是，您使用 data 变量进行 geosphere::destPoint 计算。
这是我认为应该有效的该方法的一个版本：
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
)
)
}
下一个问题是你只用卡特里娜飓风的例子定义了 1 个 x 坐标。但是，秤不知道您的半径参数，因此它们不会调整限制以适合您的半径。您可以通过设置 xmin、xmax、ymin 来规避此问题和 ymax 边界框参数，以便 scale_x_continuous() 可以了解您的半径。 （y 尺度也是如此）。我会通过对你的 ggproto 对象使用 setup_data 方法来解决这个问题。
这是我用来测试的设置数据方法，但我不是空间天才所以你必须检查这是否有意义。
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])
)
}
实施这些更改后，我得到如下图：
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
我首先定义了 class 然后是实际的 geom 函数，然而，我的输出图变得如此小型化，所以如果你能告诉我我可能哪里出错了，我将不胜感激体重秤。
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
)
)
}
)
和实际几何函数定义：
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, ...)
)
}
所以我继续绘制以下内容：
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"))
用于此目的的数据可在此处找到，作为大西洋盆地数据集 1988 - 2018： http://rammb.cira.colostate.edu/research/tropical_cyclones/tc_extended_best_track_dataset/
为了您的考虑，我使用了以下代码来整理数据：
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"))
好的，我发现了 2 个问题。问题 1 是，在您的 draw_group() ggproto 方法中，您将半径从海里转换为米（我认为），但是您将其写入 coords 变量。但是，您使用 data 变量进行 geosphere::destPoint 计算。
这是我认为应该有效的该方法的一个版本：
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
)
)
}
下一个问题是你只用卡特里娜飓风的例子定义了 1 个 x 坐标。但是，秤不知道您的半径参数，因此它们不会调整限制以适合您的半径。您可以通过设置 xmin、xmax、ymin 来规避此问题和 ymax 边界框参数，以便 scale_x_continuous() 可以了解您的半径。 （y 尺度也是如此）。我会通过对你的 ggproto 对象使用 setup_data 方法来解决这个问题。
这是我用来测试的设置数据方法，但我不是空间天才所以你必须检查这是否有意义。
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])
)
}
实施这些更改后，我得到如下图：