根据位置和日期创建星图可视化

Question

背景

我正在尝试根据 R 中给定的位置和日期创建天体地图。

理想情况下，视觉效果应如下所示：

（来源）

我确实看到了这个博客，它使用了D3 天体图，并且对于创建下面的视觉效果非常有帮助。

library(sf)
library(tidyverse)


theme_nightsky <- function(base_size = 11, base_family = "") {
  
  theme_light(base_size = base_size, base_family = base_family) %+replace% 
    theme(
      # Specify axis options, remove both axis titles and ticks but leave the text in white
      axis.title = element_blank(),
      axis.ticks = element_blank(),
      axis.text = element_text(colour = "white",size=6),
      # Specify legend options, here no legend is needed
      legend.position = "none",
      # Specify background of plotting area
      panel.grid.major = element_line(color = "grey35"),  
      panel.grid.minor = element_line(color = "grey20"),  
      panel.spacing = unit(0.5, "lines"),
      panel.background = element_rect(fill = "black", color  =  NA),  
      panel.border = element_blank(),  
      # Specify plot options
      plot.background = element_rect( fill = "black",color = "black"),  
      plot.title = element_text(size = base_size*1.2, color = "white"),
      plot.margin = unit(rep(1, 4), "lines")
    )
  
}



# Constellations Data
url1 <- "https://raw.githubusercontent.com/ofrohn/d3-celestial/master/data/constellations.lines.json"
# Read in the constellation lines data using the st_read function
constellation_lines_sf <- st_read(url1,stringsAsFactors = FALSE) %>%
                          st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>% 
                          st_transform(crs = "+proj=moll")

# Stars Data
url2 <- "https://raw.githubusercontent.com/ofrohn/d3-celestial/master/data/stars.6.json"
# Read in the stars way data using the st_read function
stars_sf <- st_read(url2,stringsAsFactors = FALSE) %>% 
            st_transform(crs = "+proj=moll")

ggplot()+
  geom_sf(data=stars_sf, alpha=0.5,color="white")+
  geom_sf(data=constellation_lines_sf, size= 1, color="white")+
  theme_nightsky()

我的问题

我设法在 R 中创建的视觉效果是（据我所知）整个天体图。我怎样才能获得我的相对位置的天体图的子集？

Answer 1

这看起来像（裁剪后的）兰伯特方位角等积投影，并且该地图似乎仅考虑了纬度（因为星图上的中心线看起来像 0 度经度）。下面的 crs 看起来差不多是正确的：

library(sf)
library(tidyverse)

toronto <- "+proj=laea +x_0=0 +y_0=0 +lon_0=0 +lat_0=43.6532"

我们需要一种方法来翻转经度坐标，使我们看起来像是在仰视天球而不是俯视天球。通过使用矩阵执行仿射变换，这相当容易做到。我们将在这里定义它：

flip <- matrix(c(-1, 0, 0, 1), 2, 2)

现在我们还需要一种方法来仅获取中心点任意方向 90 度范围内的星星（即裁剪结果）。为此，我们可以在中心点周围使用一个等于地球周长 1/4 的大缓冲区。只有与这个半球相交的星星应该是可见的。

hemisphere <- st_sfc(st_point(c(0, 43.6532)), crs = 4326) |>
              st_buffer(dist = 1e7) |>
              st_transform(crs = toronto)

我们现在可以像这样得到我们的星座：

constellation_lines_sf <- st_read(url1, stringsAsFactors = FALSE) %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>% 
  st_transform(crs = toronto) %>%
  st_intersection(hemisphere) %>%
  filter(!is.na(st_is_valid(.))) %>%
  mutate(geometry = geometry * flip) 

st_crs(constellation_lines_sf) <- toronto

而我们的明星是这样的：

stars_sf <- st_read(url2,stringsAsFactors = FALSE) %>% 
  st_transform(crs = toronto) %>%
  st_intersection(hemisphere) %>%
  mutate(geometry = geometry * flip) 

st_crs(stars_sf) <- toronto

我们还需要一个圆形蒙版来围绕最终图像进行绘制，以便生成的网格线不会延伸到圆之外：

library(grid)

mask <- polygonGrob(x = c(1, 1, 0, 0, 1, 1, 
                          0.5 + 0.46 * cos(seq(0, 2 *pi, len = 100))),
                    y =  c(0.5, 0, 0, 1, 1, 0.5, 
                           0.5 + 0.46 * sin(seq(0, 2*pi, len = 100))),
                    gp = gpar(fill = '#191d29', col = '#191d29'))

对于情节本身，我定义了一个看起来更像所需星图的主题。我将恒星星等的指数映射到大小和阿尔法，使其更加真实。

p <- ggplot() +
  geom_sf(data = stars_sf, aes(size = -exp(mag), alpha = -exp(mag)),
          color = "white")+
  geom_sf(data = constellation_lines_sf, linwidth = 1, color = "white",
          size = 2) +
  annotation_custom(circleGrob(r = 0.46, 
                               gp = gpar(col = "white", lwd = 10, fill = NA))) +
  scale_y_continuous(breaks = seq(0, 90, 15)) +
  scale_size_continuous(range = c(0, 2)) +
  annotation_custom(mask) +
  labs(caption = 'STAR MAP\nTORONTO, ON, CANADA\n9th January 2023') +
  theme_void() +
  theme(legend.position = "none",
        panel.grid.major = element_line(color = "grey35", linewidth = 1),  
        panel.grid.minor = element_line(color = "grey20", linewidth = 1),  
        panel.border = element_blank(),  
        plot.background = element_rect(fill = "#191d29", color = "#191d29"),
        plot.margin = margin(20, 20, 20, 20),
        plot.caption = element_text(color = 'white', hjust = 0.5, 
                                    face = 2, size = 25, 
                                    margin = margin(150, 20, 20, 20)))

现在，如果我们保存这个结果，请说：

ggsave('toronto.png', plot = p, width = unit(10, 'in'), 
       height = unit(15, 'in'))

我们得到