基于ggplot绘制世界地图昼夜分界线及全年动态模拟技术问询

Great question! Let's break down how to solve this step by step—combining gcIntermediate with the sunrise equation to get accurate day-night lines for any date, then building a smooth 24-hour animation.

Step 1: Understand the Sunrise Equation's Role #

The sunrise equation lets us calculate two critical values for any date/time:

• Solar declination: The angle of the Sun relative to the Earth's equator (changes with the seasons, e.g., 0° at equinoxes, ±23.44° at solstices)
• This declination defines the orientation of the day-night boundary, which is a great circle on the globe—perfect for gcIntermediate to render accurately.

Step 2: Implement the Sunrise Equation for Solar Declination #

First, write a function to compute solar declination for a given datetime:

library(lubridate)

calc_solar_declination <- function(date_time) {
  # Convert datetime to decimal Julian date
  jd <- decimal_date(date_time)
  n <- jd - 2451545.0  # Days since J2000 epoch
  
  # Calculate solar geometric parameters
  L <- (280.46646 + 0.9856473662863175 * n) %% 360  # Mean longitude
  g <- (357.52911 + 0.9856002863175 * n) %% 360      # Mean anomaly
  lambda <- L + 1.914602 * sin(g * pi/180) + 0.019993 * sin(2 * g * pi/180)  # True longitude
  
  # Compute solar declination (accounts for Earth's axial tilt)
  delta <- asin(sin(lambda * pi/180) * sin(23.44 * pi/180)) * 180/pi
  return(delta)
}

Step 3: Generate Accurate Day-Night Boundaries #

With declination in hand, we can create the great circle for the day-night line, then extend it into a closed polygon for the night shadow:

library(geosphere)

generate_daynight_shadow <- function(date_time) {
  delta <- calc_solar_declination(date_time)
  
  # Define two antipodal points on the day-night great circle
  point_north <- c(0, 90 - delta)
  point_south <- c(180, -(90 - delta))
  
  # Generate smooth great circle points
  daynight_line <- data.frame(gcIntermediate(point_north, point_south, n=200, addStartEnd=TRUE))
  colnames(daynight_line) <- c("lon", "lat")
  
  # Close the polygon to fill the night area (extend to poles)
  shadow_polygon <- rbind(
    daynight_line,
    data.frame(lon=180, lat=-90),
    data.frame(lon=-180, lat=-90),
    daynight_line[1, ]
  )
  
  shadow_polygon$datetime <- date_time
  return(shadow_polygon)
}

Step 4: Build the 24-Hour Animation #

Let's use March 21st (equinox, when declination is 0°) as an example, then generate frames for every hour and animate:

library(ggplot2)
library(ggthemes)
library(gganimate)

# Create a 24-hour UTC time sequence for March 21, 2024
start_time <- ymd_hms("2024-03-21 00:00:00", tz="UTC")
hourly_times <- start_time + hours(0:23)

# Generate shadow data for all time points
animation_data <- do.call(rbind, lapply(hourly_times, generate_daynight_shadow))

# Build the base plot
base_plot <- ggplot(animation_data) +
  borders("world", colour = "gray95", fill = "gray90") +
  geom_polygon(aes(x=lon, y=lat, group=datetime), fill="black", alpha=0.3) +
  theme_map() +
  labs(title = "Day-Night Map: {frame_time} UTC")

# Animate the plot
animate(base_plot, 
        fps=10, 
        duration=24, 
        width=800, 
        height=400,
        renderer=gifski_renderer())

Key Notes for Accuracy & Customization #

• Seasonal Adjustments: The function automatically handles different dates—try changing the start_time to June 21 (solstice) to see the shifted day-night line.
• Smoothness: Increase the n parameter in gcIntermediate for a more polished line (200 is a good balance of speed and smoothness).
• Time Zones: Stick to UTC for simplicity, but you can adjust the tz parameter if you need local time displays.
• Precision: For even more accuracy, you can add atmospheric refraction corrections to the sunrise equation, but this isn't necessary for basic day-night mapping.

内容的提问来源于stack exchange，提问作者jogall