Я работаю над созданием небесной карты на основе заданного местоположения и даты в 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, представляет собой (насколько мне известно) всю карту звездного неба. Как мне получить подмножество этой карты звездного неба для моего относительного положения?
Это выглядит как (обрезанная) азимутальная равновеликая проекция Ламберта, и карта, кажется, учитывает только широту (поскольку центральная линия выглядит как 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'))
Мы получаем
Святой дым! Отлично! Мне нужно внимательно прочитать это, чтобы все понять, но, похоже, это именно то, что мне нужно. Спасибо!
Мне просто интересно. Как мы учитываем дату и время? Это настроено в CRS со значением X или значением Y?
@Bensstats нет - небесная сфера вращается один раз в день, поэтому звездная карта будет правильной на этой широте в какой-то момент в течение 24 часов. Вы можете имитировать эффект вращения в течение 24-часового периода, добавляя или вычитая долготу из crs и центральной точки.
Хм странно. С тех пор, как я зашел на этот сайт (mapsformoments.co.uk/9-star-map?step=moment ) и сделал карту на 18 октября в Торонто, и, похоже, есть некоторая ротация (ссылка: картон. co/noy0uGrDdJx1.png). Любая подсказка, почему?
Сфера также вращается один раз в год, так что, я думаю, с тех пор произошло вращение примерно на 90 градусов? Повторяет ли изменение долготы примерно на -90 градусов октябрьский график?
Да, я только что попробовал, что там @Bensstats - если мы вычисляем пропорцию года с 18-го октября, то у нас есть -(9 + 31 + 30 + 13)/365
. Если мы затем умножим это на 360, мы получим -81,86. Если мы сделаем это нашим значением долготы в crs и центре полушария, то мы воспроизведем вашу карту от 18 октября.
Я надеюсь, что другие не столкнутся с такой же блокировкой безопасности при попытке доступа к содержимому этой ссылки на D3 Celestial Maps.