SOCR ≫ DSPA ≫ DSPA2 Topics ≫

Expanding on other DSPA2 appendices), this DSPA2 appendix covers a discussion of building a mathematical model the length/duration of daylight and demonstrating supporting data science analytics for predicting the exact length of daylight (sunset to sunrise) for a given day, longitude, latitude, and altitude.

Note: Part 2 of this DSPA Appendix 13 (Mathematical Modeling and Data Science Analytics of the Duration of Daylight) is intentionally separated as it includes large data and JavaScript code supporting the plotly app. It includes the interactive 3D temporal animation app that allows manual navigation across longitude, latitude, altitude and time (daily).

A model for computing the exact length of daylight (difference between the sunset and sunrise times) for a given day, latitude, longitude, and altitude, must utilize astronomical calculations. This requires converting the geographic position into solar coordinates and using the solar declination, \(\delta\), for a specific spacetime location and calculating the proper sunrise and sunset times.

The solar declination angle, \(\delta\), represents the latitude at which the Sun is directly overhead at solar noon. It varies throughout the year due to Earth’s tilt and orbit around the Sun. The declination can be approximated by \[\delta = 23.44^\circ \times \sin\left(\frac{2\pi}{365} \times (N - 81)\right),\] where \(N\) is the day of the year (\(N=1\) for January 1, \(N=365\) for December 31, for non leap years).

The hour angle (\(H\)) is the angle that represents the time difference from solar noon to the point at which the Sun appears on the horizon (sunrise or sunset). At sunrise and sunset, the hour angle \(H_0\) can be calculated by \[H_0 = \cos^{-1}\left(\frac{-\sin(\phi) \cdot \sin(\delta)}{\cos(\phi) \cdot \cos(\delta)}\right),\] where \(\phi\) is the latitude of the observer and \(\delta\) is the solar declination calculated earlier.

This formula assumes observation at sea level. However, at elevations above sea level, we need to adjust for altitude. For a location at altitude \(h\) (in meters), the horizon drops by an angle \(\Delta h \approx -\frac{\sqrt{2h/R}}{180/\pi}\) radians, where \(R \approx 6,371,000 \, \text{m}\) (Earth’s radius). The corrected hour angle \(H_{\text{alt}}\) becomes \[H_{\text{alt}} = \cos^{-1}\left(\frac{-\sin(\phi) \cdot \sin(\delta) - \sin(\Delta h)}{\cos(\phi) \cdot \cos(\delta)}\right) .\]

The daylight duration \(D\) in hours is twice the time from solar noon to sunset (or sunrise to solar noon), since both are symmetric around solar noon \(D = \frac{2 H_{\text{alt}}}{15} ,\) where we divide by \(15\) because each hour angle degree corresponds to \(4\) minutes.

Putting it all together, the daylight duration function \(D(\phi, N, h)\) in hours for a given latitude \(\phi\), day \(N\), and altitude \(h\) is \[D(\phi, N, h) = \frac{2}{15} \cdot \cos^{-1}\left(\frac{-\sin(\phi) \cdot \sin\left(23.44^\circ \times \sin\left(\frac{2\pi}{365} \times (N - 81)\right)\right) - \sin\left(-\frac{\sqrt{2h/R}}{180/\pi}\right)}{\cos(\phi) \cdot \cos\left(23.44^\circ \times \sin\left(\frac{2\pi}{365} \times (N - 81)\right)\right)}\right) ,\] where \(h\) is in meters (altitude above sea level) and \(R = 6,371,000\) meters is the Earth’s radius.

1 Example Calculation

For an example location at \(\phi = 40^\circ\), on the summer solstice (\(N \approx 172\)), and an altitude of \(h = 1000\) meters

  1. Compute \(\delta\).
  2. Compute \(H_{\text{alt}}\) using the altitude correction.
  3. Calculate \(D\) using the formula above.

2 R Implementation

Let’s write an R script that generates an animated 3D geomap of the earth where mouse over a geographic location (longitude, latitude, altitude) reports the daylight length for a given date that is provided by the user with a slider day between \(1\) and \(365\) (January 1 through December 31). This 3D geomap app requires combining plotly with dplyr and purrr for data processing and animation in R. The approach involves (1) Generate a grid of latitude and longitude values; (2) Calculate the daylight length for each latitude, longitude, and day of the year using the daylight duration formula; and (3) Visualize the daylight length as an animated 3D plot using plotly.

# Load necessary libraries
library(plotly)
library(dplyr)
library(purrr)

# Define the function to compute daylight duration with checks
daylight_duration <- function(latitude, day_of_year, altitude = 0) {
  # Constants
  R <- 6371000  # Earth's radius in meters
  
  # Solar declination angle
  declination <- 23.44 * sin((2 * pi / 365) * (day_of_year - 81))
  
  # Adjust for altitude in radians
  altitude_correction <- -sqrt(2 * altitude / R) * (180 / pi)
  
  # Compute the hour angle argument safely within [-1, 1]
  hour_angle_arg <- -tan(latitude * pi / 180) * tan(declination * pi / 180) - sin(altitude_correction * pi / 180)
  hour_angle_arg <- pmax(pmin(hour_angle_arg, 1), -1)
  
  # Calculate hour angle
  hour_angle <- acos(hour_angle_arg)
  
  # Calculate daylight duration in hours
  daylight_hours <- (2 * hour_angle) / (pi / 12)  # Convert radians to hours
  
  # Handle cases of 24-hour daylight or darkness near the poles
  if (hour_angle_arg <= -1) {
    daylight_hours <- 24  # Continuous daylight
  } else if (hour_angle_arg >= 1) {
    daylight_hours <- 0   # Continuous night
  }
  
  return(daylight_hours)
}

# Define the grid of latitudes and longitudes
latitude_grid <- seq(-89, 89, by = 10)
longitude_grid <- seq(-180, 180, by = 10)
days <- 1:365

# Calculate daylight duration for each (latitude, longitude, day)
data <- expand.grid(latitude=latitude_grid, longitude=longitude_grid, day=days) %>%
  mutate(daylight = map2_dbl(latitude, day, ~daylight_duration(.x, .y)))

# Create plotly animation
fig <- plot_ly(data = data, x = ~longitude, y = ~latitude, z = ~daylight,
               color = ~daylight, colors = "Blues",
               frame = ~day, type = "scatter3d", mode = "markers",
               marker = list(size = 3)) %>%
  layout(scene = list(
    xaxis = list(title = "Longitude"),
    yaxis = list(title = "Latitude"),
    zaxis = list(title = "Daylight Hours")),
    title = "Daylight Length Over Earth Throughout the Year"
    # , sliders = list(
    #   list(pad = list(t = 50),
    #     currentvalue = list(visible = TRUE, prefix = "Day of Year: "),
    #     steps = lapply(1:365, function(i) list(label = i, method = "animate", args = list(list(frame = i))))
    #   )
    # )
  )
fig # Show the plot