问题描述
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
我首先定义了类，然后定义了实际的 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,b = 90:180,]$r_se),]$alpha
)
# Creating the southwest quadrants
data_sw <- data.frame(colour = data[i,b = 180:270,]$r_sw),]$alpha
)
# Creating the northwest quadrants
data_nw <- data.frame(colour = data[i,b = 270:360,]$r_nw),]$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函数定义：
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)","yellow"))
可在此处找到用于此目的的数据，即 1988 年至 2018 年的大西洋盆地数据集： http://rammb.cira.colostate.edu/research/tropical_cyclones/tc_extended_best_track_dataset/
为了您的考虑，我使用了以下代码来整理数据：
ext_tracks_widths <- c(7,10,2,3,5,6,4,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",paste("radius_64","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"))
解决方法
好的，我发现了两个问题。问题 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,b = 90:180,d = data[i,]$r_se),]$alpha
)
# Creating the southwest quadrants
data_sw <- data.frame(colour = data[i,b = 180:270,]$r_sw),]$alpha
)
# Creating the northwest quadrants
data_nw <- data.frame(colour = data[i,b = 270:360,]$r_nw),]$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[,ymin = min(extend[,2]),ymax = max(extend[,2])
)
}
实施这些更改后，我得到了这样的数字：
