Background

I am working on trying to create a Celestial map based on a given location and date in R.

Ideally the visual would look like this:

(

I did see

My Question

The visual that I managed to create in R is (to my knowledge) the entire celestial map. How would I be able to get a subset of this celestial map for my relative position?

CodePudding user response：

This looks like a (cropped) Lambert Azimuthal equal-area projection, and the map appears to only account for latitude (since the central line looks like 0 degrees longitude on the star map). The following crs looks about right:

library(sf)
library(tidyverse)

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

We need a way to flip the longitude co-ordinates to make it appear that we are looking up at a celestial sphere rather than down on one. This is fairly easy to do by using a matrix to perform an affine transformation. We'll define that here:

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

Now we also need a way of obtaining only the stars within 90 degrees in any direction of our centre point (i.e. cropping the result). For this, we can use a large buffer around the centre point equal to 1/4 of the Earth's circumference. Only the stars that intersect this hemisphere should be visible.

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

We can now get our constellations like so:

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

And our stars like this:

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

st_crs(stars_sf) <- toronto

We'll also need a circular mask to draw around our final image, so that the resulting grid lines do not extend outside the circle:

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'))

For the plot itself, I have defined a theme that looks a bit more like the desired star map. I have mapped the exponent of star magnitude to size and alpha to make it a bit more realistic.

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)))

Now if we save this result, say with:

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

We get