dieghernan

Posted on Nov 12 Originally published at dieghernan.github.io on Jan 24, 2023

Star map with R

#rstats #astronomy #celestial #geojson

A couple of weeks ago I was doing my daily check on Stack Overflow when I found a question by Benjamin Smith that blew my mind: Creating Star Map Visualizations Based on Location and Date !!

Oh my! Can we do this with R? The answer is: Of course! In fact, Kim Fitter already worked on this two years earlier, see his Celestial Maps. So I decided to put some work into this.

Since then, I have also learned that Benjamin and I have been working in parallel on the same topic. He is preparing an R package named starBliss and hopefully this post will be of some help.

I have very little knowledge of this topic, so if you find errors or have any suggestions, please let me know in the comments section below.

First things first: the data

The initial data source for all these projects (Kim's, Benjamin's and mine) is the same, and it is provided by the D3 plugin d3-celestial by Olaf Frohn. As Kim Fitter pointed out in his post, these data files present the problem (experienced by almost any sf user) of lines crossing the international date line (longitude 180°). I also found that some files are not valid as per sf::st_make_valid().

Solution? I processed and fixed almost every file, including the files corresponding to the Milky Way and the lines for Chinese constellations, manually, to provide a set of files. That is the origin of my project Celestial Data, which provides all these files in several spatial formats. Please check out the repo to learn more about it.

Creating a star map with R

The first step is to load libraries that will help us with this cosmic task:

# Spatial manipulation
library(sf)
library(s2)
library(nominatimlite)

## Wrangling data and dates
library(dplyr)
library(lubridate)
library(lutz)

## Visualization
library(ggplot2)
library(ggfx)
library(ggshadow)

Helper functions

Now we prepare some helper functions:

load_celestial() just downloads the corresponding .geojson from the Celestial Data repo¹ to a specific directory cachedir and loads it with sf::st_read().
pretty_lonlat() is a labeler that returns a decimal longitude or latitude coordinate in the format degrees/minutes/seconds (e.g., a latitude of 34.72782 would be converted into 34° 43' 40.15" N).

Show load_celestial and pretty_lonlat()

load_celestial <- function(filename,
                           url = "https://cdn.jsdelivr.net/gh/dieghernan/celestial_data@main/data/",
                           cachedir = tempdir()) {
  if (!dir.exists(cachedir)) {
    stop(
      "Please create ",
      path.expand(cachedir),
      " directory",
      "first"
    )
  }

  url <- file.path(url, filename)
  local_path <- file.path(cachedir, filename)


  if (!file.exists(local_path)) {
    download.file(url, local_path, mode = "wb", quiet = TRUE)
  }

  celestial <- sf::st_read(local_path, quiet = TRUE)

  return(celestial)
}

pretty_lonlat <- function(x, type, accuracy = 2) {
  positive <- x >= 0

  # Decompose
  x <- abs(x)
  D <- as.integer(x)
  m <- (x - D) * 60
  M <- as.integer(m)
  S <- round((m - M) * 60, accuracy)

  # Get label
  if (type == "lon") {
    lab <- ifelse(positive > 0, "E", "W")
  } else {
    lab <- ifelse(positive > 0, "N", "S")
  }


  # Compose
  label <- paste0(D, "\u00b0 ", M, "' ", S, '\" ', lab)
  return(label)
}

Additionally, you may notice that on this d3-celestial demo there is some degree of rotation depending on the location and the time. I found how this is done on the d3-celestial plugin and I found the function getMST(dt, lng), which I ported to R (get_mst()). The function uses Meeus's (1998) formula for Greenwich mean sidereal time, adds the observer's east-positive longitude and wraps the result to the interval [0°, 360°). This gives local mean sidereal time in degrees. For more detail, see this post by James Still.

The inputs are a POSIXct date-time and the observer's longitude in degrees. The example uses the returned angle to rotate the star map's projection.

Show get_mst()

# Derive rotation degrees of the projection given a date and a longitude
get_mst <- function(dt, lng) {
  desired_date_utc <- lubridate::with_tz(dt, "UTC")


  yr <- lubridate::year(desired_date_utc)
  mo <- lubridate::month(desired_date_utc)
  dy <- lubridate::day(desired_date_utc)
  h <- lubridate::hour(desired_date_utc)
  m <- lubridate::minute(desired_date_utc)
  s <- lubridate::second(desired_date_utc)

  if ((mo == 1) || (mo == 2)) {
    yr <- yr - 1
    mo <- mo + 12
  }

  # Adjust times before Gregorian Calendar
  # See https://squarewidget.com/julian-day/
  if (lubridate::as_date(dt) > as.Date("1582-10-14")) {
    a <- floor(yr / 100)
    b <- 2 - a + floor(a / 4)
  } else {
    b <- 0
  }
  c <- floor(365.25 * yr)
  d <- floor(30.6001 * (mo + 1))

  # days since J2000.0
  jd <- b + c + d - 730550.5 + dy + (h + m / 60 + s / 3600) / 24
  jt <- jd / 36525

  # Rotation
  mst <- 280.46061837 + 360.98564736629 * jd +
    0.000387933 * jt^2 - jt^3 / 38710000.0 + lng

  # Modulo 360 degrees
  mst <- mst %% 360

  return(mst)
}

The final result would have a spherical outline. That means that we would need to perform a spherical cut. Did you know that in r-spatial the Earth is no longer flat? Thanks to s2 we can overcome this issue. Additionally, we would get rid of artifacts derived from the changes on the projection. This also includes some refinements to avoid empty/non-valid geometries as well as GEOMETRYCOLLECTION handling. The function sf_spherical_cut() would do that for us.

Show sf_spherical_cut()

# Cut a sf object with a buffer using spherical s2 geoms
# Optionally, project and flip

sf_spherical_cut <- function(x, the_buff, the_crs = sf::st_crs(x), flip = NULL) {
  # Get geometry type
  geomtype <- unique(gsub("MULTI", "", sf::st_geometry_type(x)))[1]

  # Keep the data frame, s2 drops it
  the_df <- sf::st_drop_geometry(x)
  the_geom <- sf::st_geometry(x)
  # Convert to s2 if needed
  if (!inherits(the_buff, "s2_geography")) {
    the_buff <- sf::st_as_s2(the_buff)
  }

  the_cut <- the_geom %>%
    # Cut with s2
    sf::st_as_s2() %>%
    s2::s2_intersection(the_buff) %>%
    # Back to sf and add the df
    sf::st_as_sfc() %>%
    sf::st_sf(the_df, geometry = .) %>%
    dplyr::filter(!sf::st_is_empty(.)) %>%
    sf::st_transform(crs = the_crs)

  # If it is not POINT filter by valid and non-empty
  # This is for performance
  if (!geomtype == "POINT") {
    # If any is GEOMETRYCOLLECTION extract the right value
    if (any(sf::st_geometry_type(the_cut) == "GEOMETRYCOLLECTION")) {
      the_cut <- the_cut %>%
        sf::st_collection_extract(type = geomtype, warn = FALSE)
    }

    the_cut <- the_cut %>%
      dplyr::filter(!is.na(sf::st_is_valid(.)))
  }

  if (!is.null(flip)) {
    the_cut <- the_cut %>%
      dplyr::mutate(geometry = geometry * flip) %>%
      sf::st_set_crs(the_crs)
  }

  return(the_cut)
}

Inputs

Now we are ready to start creating our visualization. We need only two inputs:

A desired location, which we will geocode with nominatimlite.
A specific moment in time.

# Inputs
desired_place <- "Madrid, Spain"

# We are not using yet the timezone
desired_date <- make_datetime(
  year = 2015,
  month = 9,
  day = 22,
  hour = 3,
  min = 45
)

# Geocode place with nominatimlite
desired_place_geo <- geo_lite(desired_place, full_results = TRUE)

desired_place_geo %>%
  select(address, lat, lon)
#> # A tibble: 1 × 3
#>   address                                                                                    lat   lon
#>   <chr>                                                                                    <dbl> <dbl>
#> 1 Madrid, Área metropolitana de Madrid y Corredor del Henares, Comunidad de Madrid, España  40.4 -3.70

# And get the coordinates
desired_loc <- desired_place_geo %>%
  select(lat, lon) %>%
  unlist()

desired_loc
#>       lat       lon
#> 40.416705 -3.703582

With respect to our object desired_date, it is quite relevant for accurate plotting to specify the correct time zone. Since we already know the latitude and longitude of our desired location, we can easily get that with the lutz package:

desired_date
#> [1] "2015-09-22 03:45:00 UTC"

# Get tz
get_tz <- tz_lookup_coords(desired_loc[1], desired_loc[2], warn = FALSE)

get_tz
#> [1] "Europe/Madrid"

# Force it to be local time
desired_date_tz <- force_tz(desired_date, get_tz)

desired_date_tz
#> [1] "2015-09-22 03:45:00 CEST"

About time zones

Some online shops that create this kind of map (I won't post links) include this script:

...

'selectedHour': '22',
'selectedMinute': '00',

...

This means that those shops are really creating the map at YYYY-MM-DD 22:00:00 UTC. If you want to replicate that exactly (even though the resulting night sky is not accurate—for example, in New Zealand the local time at that moment would be 10:00, so no stars would be visible), you would need to adjust desired_date_tz as follows:

as_datetime(paste(as.Date(desired_date_tz), "22:00:00"), tz = "UTC")
#> [1] "2015-09-22 22:00:00 UTC"

# That would really correspond to 10:00
as_datetime(paste(as.Date(desired_date_tz), "22:00:00"), tz = "UTC") %>%
  with_tz("Pacific/Auckland")
#> [1] "2015-09-23 10:00:00 NZST"

Setup

Now we can start creating our buffers and projections, which will help us to crop the celestial data objects.

I also noticed that the location demo of d3-celestial.js uses Airy projection, so we are going to replicate that as well:

# Get the rotation and prepare buffer and projection

# Get right degrees
lon_prj <- get_mst(desired_date_tz, desired_loc[2])
lat_prj <- desired_loc[1]

c(lon_prj, lat_prj)
#>      lon      lat
#> 23.15892 40.41670

# Create proj4string w/ Airy projection

target_crs <- paste0("+proj=airy +x_0=0 +y_0=0 +lon_0=", lon_prj, " +lat_0=", lat_prj)


target_crs
#> [1] "+proj=airy +x_0=0 +y_0=0 +lon_0=23.1589164999314 +lat_0=40.4167047"

# Flip celestial objects to give the impression of viewing the sky from Earth,
# rather than viewing the Earth from the sky.
# https://stackoverflow.com/a/75064359/7877917
# Flip matrix for affine transformation

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


# And create an s2 buffer of the visible hemisphere at the given location
hemisphere_s2 <- s2_buffer_cells(
  as_s2_geography(
    paste0("POINT(", lon_prj, " ", lat_prj, ")")
  ),
  9800000,
  max_cells = 5000
)

# This one is for plotting
hemisphere_sf <- hemisphere_s2 %>%
  st_as_sf() %>%
  st_transform(crs = target_crs) %>%
  st_make_valid()

Celestial data

Now, we can load the data of our choice. In this case, I have chosen to represent the Milky Way, constellation lines and stars.

We also add some additional variables that would help us to improve the visualization.

mw <- load_celestial("mw.min.geojson")

# Add colors to MW to use on fill
cols <- colorRampPalette(c("white", "yellow"))(5)
mw$fill <- factor(cols, levels = cols)

ggplot(mw) +
  geom_sf(aes(fill = fill)) +
  scale_fill_identity()


# And process it

# Cut to buffer
mw_end <- sf_spherical_cut(mw,
  the_buff = hemisphere_s2,
  # Change the crs
  the_crs = target_crs,
  flip = flip_matrix
)


ggplot(mw_end) +
  geom_sf(aes(fill = fill)) +
  scale_fill_identity()

Now it is the turn of the constellations:

const <- load_celestial("constellations.lines.min.geojson")

ggplot(const) +
  geom_sf() +
  coord_sf(expand = FALSE)


# Cut to buffer

const_end <- sf_spherical_cut(const,
  the_buff = hemisphere_s2,
  # Change the crs
  the_crs = target_crs,
  flip = flip_matrix
)


ggplot(const_end) +
  geom_sf() +
  coord_sf(expand = FALSE)

And finally the stars:

stars <- load_celestial("stars.6.min.geojson")

ggplot(stars) +
  # We use relative brightness (br) as aes
  geom_sf(aes(size = br, alpha = br), shape = 16) +
  scale_size_continuous(range = c(0.5, 6)) +
  scale_alpha_continuous(range = c(0.1, 0.8)) +
  coord_sf(expand = FALSE)


# Cut to buffer

stars_end <- sf_spherical_cut(stars,
  the_buff = hemisphere_s2,
  # Change the crs
  the_crs = target_crs,
  flip = flip_matrix
)

ggplot(stars_end) +
  # We use relative brightness (br) as aes
  geom_sf(aes(size = br, alpha = br), shape = 16) +
  scale_size_continuous(range = c(0.5, 6)) +
  scale_alpha_continuous(range = c(0.1, 0.8))

Graticules

We are also going to include graticules, so the Earth's poles can be quickly spotted. In this case we don't apply any affine transformation, so the flip parameter of sf_spherical_cut() needs to be set as NULL.

grat <- st_graticule(
  ndiscr = 5000,
  lat = seq(-90, 90, 10),
  lon = seq(-180, 180, 30)
)

ggplot(grat) +
  geom_sf() +
  coord_sf(expand = FALSE)


# Cut to buffer; we don't flip this one (it is not an object of the space)
grat_end <- sf_spherical_cut(
  x = grat,
  the_buff = hemisphere_s2,
  # Change the crs
  the_crs = target_crs
)


ggplot(grat_end) +
  geom_sf() +
  coord_sf(expand = FALSE)

Visualization with ggplot2

We are almost set! For preparing the final map, first we are going to create the corresponding labels, which will be included as caption on the ggplot2 map:

lat_lab <- pretty_lonlat(desired_loc[1], type = "lat")
lon_lab <- pretty_lonlat(desired_loc[2], type = "lon")

pretty_labs <- paste(lat_lab, "/", lon_lab)

cat(pretty_labs)
#> 40° 25' 0.14" N / 3° 42' 12.9" W

# Create final caption to put on bottom

pretty_time <- paste(
  # Pretty Day
  scales::label_date(
    format = "%d %b %Y",
    locale = "en"
  )(desired_date_tz),
  # Pretty Hour
  format(desired_date_tz, format = "%H:%M", usetz = TRUE)
)

cat(pretty_time)
#> 22 Sep 2015 03:45 CEST

# Our final caption
caption <- toupper(paste0(
  "Star Map\n",
  desired_place, "\n",
  pretty_time, "\n",
  pretty_labs
))


cat(caption)
#> STAR MAP
#> MADRID, SPAIN
#> 22 SEP 2015 03:45 CEST
#> 40° 25' 0.14" N / 3° 42' 12.9" W

We can enhance the visualization by applying some interesting effects:

We want the Milky Way to appear a bit blurry instead of as a Well-Known geometry. With this effect we can mimic how we really see it from the Earth. So we can use ggfx::with_blur() to get this effect.
We can also add a glowing effect to our stars and constellations (have a look at Dominic Royé's post Firefly Cartography to know more about this). The only drawback is that I was not able to use ggshadow with LINESTRING (Dominic shows how to do it with POINT), so instead I converted my lines (the constellations) to coordinates and applied ggshadow::geom_glowpath().

So we are ready now to create the final visualization:

# Prepare MULTILINESTRING

const_end_lines <- const_end %>%
  st_cast("MULTILINESTRING") %>%
  st_coordinates() %>%
  as.data.frame()


ggplot() +
  # Graticules
  geom_sf(data = grat_end, color = "grey60", linewidth = 0.25, alpha = 0.3) +
  # A blurry Milky Way
  with_blur(
    geom_sf(
      data = mw_end, aes(fill = fill), alpha = 0.1, color = NA,
      show.legend = FALSE
    ),
    sigma = 8
  ) +
  scale_fill_identity() +
  # Glowing stars
  geom_glowpoint(
    data = stars_end, aes(
      alpha = br, size =
        br, geometry = geometry
    ),
    color = "white", show.legend = FALSE, stat = "sf_coordinates"
  ) +
  scale_size_continuous(range = c(0.05, 0.75)) +
  scale_alpha_continuous(range = c(0.1, 0.5)) +
  # Glowing constellations
  geom_glowpath(
    data = const_end_lines, aes(X, Y, group = interaction(L1, L2)),
    color = "white", size = 0.5, alpha = 0.8, shadowsize = 0.4, shadowalpha = 0.01,
    shadowcolor = "white", linejoin = "round", lineend = "round"
  ) +
  # Border of the sphere
  geom_sf(data = hemisphere_sf, fill = NA, color = "white", linewidth = 1.25) +
  # Caption
  labs(caption = caption) +
  # And end with theming
  theme_void() +
  theme(
    text = element_text(colour = "white"),
    panel.border = element_blank(),
    plot.background = element_rect(fill = "#191d29", color = "#191d29"),
    plot.margin = margin(20, 20, 20, 20),
    plot.caption = element_text(
      hjust = 0.5, face = "bold",
      size = rel(1),
      lineheight = rel(1.2),
      margin = margin(t = 40, b = 20)
    )
  )

Voilà! I checked the results several times against those provided by d3-celestial.js on the location demo and the underlying calculations in JavaScript, and everything seems to be up and running.

Extra: Chinese constellations

Celestial Data also provides data for traditional Chinese constellations, so we can create a similar map with this entirely different set of geometries:

const_cn <- load_celestial("constellations.lines.cn.min.geojson")

# Cut and prepare for geom_glowpath() on a single step
const_cn_end_lines <- sf_spherical_cut(const_cn,
  the_buff = hemisphere_s2,
  # Change the crs
  the_crs = target_crs,
  flip = flip_matrix
) %>%
  # To paths
  st_cast("MULTILINESTRING") %>%
  st_coordinates() %>%
  as.data.frame()


ggplot() +
  # Graticules
  geom_sf(data = grat_end, color = "grey60", linewidth = 0.25, alpha = 0.3) +
  # A blurry Milky Way
  with_blur(
    geom_sf(
      data = mw_end, aes(fill = fill), alpha = 0.1, color = NA,
      show.legend = FALSE
    ),
    sigma = 8
  ) +
  scale_fill_identity() +
  # Glowing stars
  geom_glowpoint(
    data = stars_end, aes(
      alpha = br, size =
        br, geometry = geometry
    ),
    color = "white", show.legend = FALSE, stat = "sf_coordinates"
  ) +
  scale_size_continuous(range = c(0.05, 0.75)) +
  scale_alpha_continuous(range = c(0.1, 0.5)) +
  # Glowing constellations
  geom_glowpath(
    data = const_cn_end_lines, aes(X, Y, group = interaction(L1, L2)),
    color = "white", size = 0.5, alpha = 0.8, shadowsize = 0.4, shadowalpha = 0.01,
    shadowcolor = "white", linejoin = "round", lineend = "round"
  ) +
  # Border of the sphere
  geom_sf(data = hemisphere_sf, fill = NA, color = "white", linewidth = 1.25) +
  # Caption
  labs(caption = caption) +
  # And end with theming
  theme_void() +
  theme(
    text = element_text(colour = "white"),
    panel.border = element_blank(),
    plot.background = element_rect(fill = "#191d29", color = "#191d29"),
    plot.margin = margin(20, 20, 20, 20),
    plot.caption = element_text(
      hjust = 0.5, face = "bold",
      size = rel(1),
      lineheight = rel(1.2),
      margin = margin(t = 40, b = 20)
    )
  )

References

Fitter, K. (2019, February 21). Celestial maps. Kim Fitter. https://kimnewzealand.github.io/2019/02/21/celestial-maps/
Frohn, O. (2015). d3-celestial [Computer software]. GitHub. https://github.com/ofrohn/d3-celestial/
Frohn, O., & Hernangómez, D. (2023). Celestial data (Version 1.0.1) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.7561601 (Project documentation)
Meeus, J. (1998). Astronomical algorithms (2nd ed.). Willmann-Bell.
Pebesma, E., & Dunnington, D. (2020, June 17). In r-spatial, the Earth is no longer flat. r-spatial. https://r-spatial.org/r/2020/06/17/s2.html
Royé, D. (2021, June 1). Firefly cartography. Dr Dominic Royé. https://dominicroye.github.io/blog/firefly-maps/
Still, J. (2020, March 22). Astronomical calculations: Sidereal time. James Still. https://squarewidget.com/astronomical-calculations-sidereal-time/

For the datasets used in this tutorial, see Celestial data.

In fact, the download is performed via jsDelivr, which distributes files hosted on GitHub via CDN. This is supposed to improve performance, but in any case the underlying data source is the GitHub repo. ↩

DEV Community