Speed and Accuracy Tests
Package ‘SunCalcMeeus’ 0.1.2
Pedro J. Aphalo
2025-01-09
Source:vignettes/articles/tests.Rmd
tests.RmdIntroduction
R package ‘SunCalcMeeus’ collects functions related to the computation of the position of the sun and times, including day length and night length, from package ‘photobiology’. It also includes function for computing empirical estimates of the relative air mass (AM) for different sun elevation angles.
The position of the sun is computed using Meeus well known algorithm, which is accurate over a broad range of dates both into the past and future. The algorithm is described as useful for dates between 1000 BC to 3000 AD, however accuracy of the implementation in ‘SunCalcMeeus’ could be affected by how Julian days are computed and the limited precision of the computer representation of real numbers as double precision floating point numbers. Results for years 1700 to 2100 are checked below against the tables of the US Naval Observatory as reference.
Note: Although for many uses, including biology, energy use and history research and applications, Meeus algorithm is considered very accurate, for use in astronomy, it is considered not the most accurate. One collection of libraries for astronomy, in FORTRAN and ANSI C versions, is provided by the International Astronomy Union under the name of Standards of Fundametal Astronomy SOFA. The NASA web page JPL Horizons provides on-line computed up-to-date ephemerides for for huge number of celestial bodies covering from 9999 BC to 9999 AD.
A comparison to package ‘suncalc’ is also done both for performance and accuracy. This is of interest as the two packages are likely to differ both in performance and accuracy, and their use could preferable in different situations.
Accuracy
Sunset and sunrise times
In this section I check the values from ‘SunCalcMeeus’ against the
on-line data from the Astronomical Applications
Department of the U.S. Naval Observatory. For any comparison of
sunrise or sunset times, it is important to ensure that the twilight
definition used is the same. In other words, at what instant in time the
start and end of the day are assumed to be is crucial. In addition, if
an atmospheric refraction correction is applied or not must also be
consistent. The astronomical definition of sunrise and sunset is when
the upper rim of the solar disk is at the horizon, and most frequently a
refraction correction is applied when the observer is on Earth. This
corresponds to twilight = "sunlight" in the functions from
‘SunCalcMeeus’. Some uncertainty remains as the refraction correction is
affected by the water column in the atmosphere and temperature, the
assumed values may differ.
The format of the table downloaded, or more precisely copied and paste from the web page, is such that converting the data into a shape compatible with that used by ‘SunCalcMeeus’ takes a some effort. I use for the comparison the geographic location of Greenwhich, U.K., and the GMT time zone.
# map
vars2months <- rep(month.abb, each = 2)
names(vars2months) <- paste("X", 2L:25L, sep = "")
vars2rise_or_set <- rep(c("sunrise", "sunset"), times = 12)
names(vars2rise_or_set) <- paste("X", 2L:25L, sep = "")
for (this.year in c(1700, 1900, 2000, 2100)) {
read_fwf(paste("data/USNO-", this.year, ".txt", sep = ""), skip = 9, n_max = 31,
col_types = "ccccccccccccccccccccccccc") |>
pivot_longer(X2:X25, values_to = "time", cols_vary = "slowest") |>
mutate(year = this.year,
month = vars2months[name],
day = X1,
event = vars2rise_or_set[name],
X1 = NULL,
name = NULL,
hhmm = ymd_hm(paste(year, "-", month, "-", day, " ", time, sep = ""),
tz = "UTC"),
day = ymd(paste(year, "-", month, "-", day, sep = ""),
tz = "UTC"),
time = NULL)|>
select(day, event, hhmm) -> all.data
sunrises.df <- subset(all.data, event == "sunrise", select = -2)
names(sunrises.df)[2] <- "usno.sunrise"
sunsets.df <- subset(all.data, event == "sunset")
sunrises.df$usno.sunset <- sunsets.df$hhmm
usno.df <- na.omit(sunrises.df)
meeus_times.df <- day_night(usno.df$day, unit.out = "datetime", twilight = "sunlight") # solar rim + refraction correction
full_join(meeus_times.df, usno.df, by = join_by(day)) |>
mutate(sunrise.diff = sunrise - usno.sunrise,
sunset.diff = sunset - usno.sunset) -> temp.df
assign(paste("times_", this.year, ".df", sep = ""), temp.df)
}
times_all.df <- bind_rows(times_1700.df, times_1900.df, times_2000.df, times_2100.df)
times_all.df$year <- year(times_all.df$day)The quantiles for the pooled differences from all years tested are:
quantile(times_all.df$sunrise.diff)## Time differences in secs
## 0% 25% 50% 75% 100%
## -69.1767883 -22.6055547 0.7093608 23.3451028 58.2047347
quantile(times_all.df$sunset.diff)## Time differences in secs
## 0% 25% 50% 75% 100%
## -60.0021386 -21.5505843 -0.3213501 25.5355292 70.3363082Differences were in all cases less than s, while the resolution of the data from the U.S. Naval Observatory is 1 min. This good match was true for years 1700, 1900, 2000 and 2100, and for both sunrise and sunset times. So, Meeus algorithm is generating values differing little from known-good values.
We plot the differences for these four years, both for sunrise and sunset. First as a density distribution and farther down in a scatter plots.
ggplot(times_all.df, aes(sunrise.diff)) +
stat_density() +
facet_wrap(facets = vars(year), labeller = label_both) +
expand_limits(x = c(-80, 80)) +
labs(x = "Difference in sunrise time (s)")
ggplot(times_all.df, aes(sunset.diff)) +
stat_density() +
facet_wrap(facets = vars(year), labeller = label_both) +
expand_limits(x = c(-80, 80)) +
labs(x = "Difference in sunset time (s)")
ggplot(times_all.df, aes(day, sunrise.diff)) +
geom_point(size = 0.7, alpha = 0.5) +
stat_smooth(method = "gam") +
labs(y = "Difference in sunrise time (s)", x = "Date") +
scale_x_datetime(date_labels = "%b") +
facet_wrap(facets = vars(year), scales = "free_x", labeller = label_both)
ggplot(times_all.df, aes(day, sunset.diff)) +
geom_point(size = 0.7, alpha = 0.5) +
stat_smooth(method = "gam") +
labs(y = "Difference in sunset time (s)", x = "Date") +
scale_x_datetime(date_labels = "%b") +
facet_wrap(facets = vars(year), scales = "free_x", labeller = label_both)
The scatterplots look like bands approximately 60 s wide, this is likely due, at least in part to the rounding to the nearest whole minute in the data used as reference. Such rounding can be expected to result in deviations within s of the full-resolution values. The S-shaped underlying line is a limitation of the algorithm or its implementation.
‘SunCalcMeeus’ provides accurate estimates, but how does its performance compare to ‘suncalc’? Is the possibly slower performance worth the expected improvement in accuracy?
Performance
Meeus mathematical algorithm is more computationally intensive than other less accurate ones, such as that implemented in package ‘suncalc’. However, I have optimized the implementation in ‘SunCalcMeeus’ for the case of long vectors of instants in time or of dates at each geographic location.
Note: The equations used in ‘suncalc’ are “mostly” from the web page Astronomy Answers. I was unable to find a published source for the algorithm other than the web page. It is not clear from the documentation if an atmospheric refraction correction is applied or not. The refraction correction slightly lowers the elevation angle at sunset and sunrise, and would introduce an error that is constant from year to year.
Sun position
To check the call overhead, I start by benchmarking the computation of the sun position at a single instant in time at a single geographic location.
time <- now()
microbenchmark(sun_angles(time),
getSunlightPosition(date = time, lat = 51.5, lon = 0),
times = 500L,
unit = "millisecond")## Unit: milliseconds
## expr min lq
## sun_angles(time) 3.905001 4.042152
## getSunlightPosition(date = time, lat = 51.5, lon = 0) 0.686902 0.732951
## mean median uq max neval cld
## 4.6711798 4.152152 4.349352 13.9343 500 a
## 0.9702228 0.776201 0.869351 45.9105 500 bI next test the performance of the computation of the position of the
sun with ‘SunCalcMeeus’ using 10 000 instants in time stored in a
POSIXct vector and nine geographic locations in a
data.frame with columns for longitude and latitude and rows
for geographic locations.
times <- ymd_hms("2024-12-11 00:00:00") + minutes(1:1e4)
locations <- data.frame(lon = c(0, 0, 0, 90, 90, 90, -90, -90, -90),
lat = c(0, 30, 60, 0, 30, 60, 0, 30, 60))
range(times)## [1] "2024-12-11 00:01:00 UTC" "2024-12-17 22:40:00 UTC"We benchmark performance by evaluating the statement 20 times. The median time is the most informative.
microbenchmark(sun_angles(times, geocode = locations),
times = 20L)## Unit: milliseconds
## expr min lq mean median
## sun_angles(times, geocode = locations) 205.9283 221.1198 303.0577 304.1168
## uq max neval
## 376.31 414.4792 20Next I test the performance with 1 000 000 instants in time and one geographic location, the default one (Greewich, UK).
## [1] "2024-12-11 00:01:00 UTC" "2026-11-05 10:40:00 UTC"There is one problem, possibly a division by zero, affecting azimuth for one instant in time. I try to locate to which instant in time it corresponds.
sun_angles.df <- sun_angles(times)## Warning in acos(((sin(lat.rad) * cos(zenith.angle.rad)) -
## sin(declin.rad))/(cos(lat.rad) * : NaNs produced
## Warning in acos(((sin(lat.rad) * cos(zenith.angle.rad)) -
## sin(declin.rad))/(cos(lat.rad) * : NaNs produced
sapply(sun_angles.df, anyNA)## time tz solartime longitude latitude address
## FALSE FALSE FALSE FALSE FALSE FALSE
## azimuth elevation declination eq.of.time hour.angle distance
## TRUE FALSE FALSE FALSE FALSE FALSE## [1] 948230## # A tibble: 1 × 12
## time tz solartime longitude latitude address azimuth
## <dttm> <chr> <solar_tm> <dbl> <dbl> <chr> <dbl>
## 1 2026-09-30 11:50:00 UTC 11:59:59 0 51.5 Greenwich NaN
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## # hour.angle <dbl>, distance <dbl>I run the benchmarks with the longer vector of instants in time:
microbenchmark(sun_angles(times),
getSunlightPosition(date = times, lat = 51.5, lon = 0),
times = 5L,
unit = "second")## Unit: seconds
## expr min lq
## sun_angles(times) 2.2363986 2.2606236
## getSunlightPosition(date = times, lat = 51.5, lon = 0) 0.7485143 0.7603652
## mean median uq max neval cld
## 2.3533759 2.316995 2.3392938 2.613568 5 a
## 0.8594519 0.882598 0.9056159 1.000166 5 bsun_angles() takes between 2.5 and 3.2 times as much
time to run than getSunlightPosition() but returns six
computed values while getSunlightPosition returns only two.
For most uses they are both most likely fast enough.
Sunrise and sunset times
Checking performance similarly as above, but day by day for 600 years.
## [1] 219146
range(dates)## [1] "1700-01-01 UTC" "2300-01-01 UTC"First just a single date to asses the call overhead.
microbenchmark(day_night(dates[1]),
day_night(dates[1], unit.out = "datetime"),
getSunlightPosition(date = dates[1], lat = 51.5, lon = 0),
times = 100L,
unit = "millisecond")## Unit: milliseconds
## expr min lq
## day_night(dates[1]) 4.033501 4.198301
## day_night(dates[1], unit.out = "datetime") 5.163101 5.454550
## getSunlightPosition(date = dates[1], lat = 51.5, lon = 0) 0.714001 0.851651
## mean median uq max neval cld
## 4.952900 4.448851 4.91990 16.854700 100 a
## 6.703946 5.694251 6.70000 17.997201 100 b
## 1.100248 0.931001 1.07935 4.890401 100 cAnd with every day in 600 years.
microbenchmark(day_night(dates),
day_night(dates, unit.out = "datetime"),
getSunlightPosition(date = dates, lat = 51.5, lon = 0),
times = 10L,
unit = "second")## Unit: seconds
## expr min lq
## day_night(dates) 0.4289319 0.4551011
## day_night(dates, unit.out = "datetime") 0.6713478 0.6753242
## getSunlightPosition(date = dates, lat = 51.5, lon = 0) 0.1450995 0.1464530
## mean median uq max neval cld
## 0.5102445 0.4898529 0.5368722 0.6432271 10 a
## 0.8002286 0.7407470 0.9206099 1.0010827 10 b
## 0.1749730 0.1507766 0.1674113 0.3476669 10 cThe difference is similar to that above for the position. But once again, for most uses they are both fast enough.
meeus_dn.df <- day_night(dates[1], tz = "UTC", unit.out = "datetime", twilight = "sunlight")
colnames(meeus_dn.df)## [1] "day" "tz" "twilight.rise" "twilight.set"
## [5] "longitude" "latitude" "address" "sunrise"
## [9] "noon" "sunset" "daylength" "nightlength"Now we use the ‘suncalc’ package to compute values for the same day.
suncalc_dn.df <- getSunlightTimes(date = as.Date(dates[1], tz = "UTC"), lat = 51.5, lon = 0)
colnames(suncalc_dn.df)## [1] "date" "lat" "lon" "solarNoon"
## [5] "nadir" "sunrise" "sunset" "sunriseEnd"
## [9] "sunsetStart" "dawn" "dusk" "nauticalDawn"
## [13] "nauticalDusk" "nightEnd" "night" "goldenHourEnd"
## [17] "goldenHour"Consistency between ‘suncalc’ and ‘SunCalcMeeus’
Sun position
The next question is do the values returned by functions from both packages agree? and does the agreement hold for past and future centuries. Let’s check every hour for 600 years starting from 1 January 1700.
## [1] "1700-01-01 01:00:00 UTC" "2300-01-01 00:00:00 UTC"
meeus.df <- sun_angles(times, use.refraction = TRUE)
meeus.df## # A tibble: 5,259,480 × 12
## time tz solartime longitude latitude address azimuth
## <dttm> <chr> <solar_tm> <dbl> <dbl> <chr> <dbl>
## 1 1700-01-01 01:00:00 UTC 00:55:49 0 51.5 Greenwich 26.0
## 2 1700-01-01 02:00:00 UTC 01:55:48 0 51.5 Greenwich 49.0
## 3 1700-01-01 03:00:00 UTC 02:55:47 0 51.5 Greenwich 66.7
## # ℹ 5,259,477 more rows
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## # hour.angle <dbl>, distance <dbl>Now we use the ‘suncalc’ package to compute the same values.
suncalc.df <- getSunlightPosition(date = times, lat = 51.5, lon = 0)
suncalc.df$elevation <- suncalc.df$altitude * 180 / pi # radians to degrees
suncalc.df$azimuth <- suncalc.df$azimuth * 180 / pi # radians to degrees
suncalc.df$azimuth <- suncalc.df$azimuth + 180 # radians to degrees
head(suncalc.df, n = 3L)## date lat lon altitude azimuth elevation
## 1 1700-01-01 01:00:00 51.5 0 -1.0442259 19.96054 -59.82974
## 2 1700-01-01 02:00:00 51.5 0 -0.9572501 44.09387 -54.84639
## 3 1700-01-01 03:00:00 51.5 0 -0.8262301 62.75066 -47.33950We compute the differences in sun elevations, and their range. Over this extended period, the maximum and minimum differences are substantial. We plot a smooth spline on the absolute value of the differences, as from the range we already know that they maxima and minimum differences are of opposite signs. The sun elevation values from the two functions are closest during the first half of the 21st century. On the other hand we see that the fitted spline does not reach the extreme values by a large margin. What is going on?
differences <- suncalc.df$elevation - meeus.df$elevation
differences.range <- range(differences)
differences.range## [1] -3.672678 3.069026
ggplot(data.frame(times, differences), aes(times, abs(differences))) +
stat_smooth(method = "gam") +
stat_panel_counts(label.x = "centre") +
labs(x = "Time", y = "Absolute difference in sun elevation (degrees)") +
expand_limits(y = 0)
If we plot the 2D density distribution we get a very surprising pattern. So as we move into the past we mainly have three concentrations of differences: positive, no difference or negative, with fewer intermediate values. When we move into the future we get four rays of concentrated difference values.
ggplot(data.frame(times, differences), aes(times, differences)) +
stat_hdr() +
stat_panel_counts(label.x = "centre") +
labs(x = "Time", y = "Difference in sun elevation (degrees)")
A difference in sun elevation of 3.5 degrees is too much for many applications. So, with Meeus algorithm we obtain consistenty accurate values across 600 years, while with ‘suncalc’ the accuracy degrades as we move away from the current half century.
We can expect this drift to also affect times. In the next section I check the difference in the estimates of time related quantities.
Sunrise and sunset times
Lets check, as above, every day for 600 years starting from 1 January 1700.
## [1] "1700-01-01 UTC" "2300-01-01 UTC"
meeus_dn.df <- day_night(times, tz = "UTC", unit.out = "datetime", twilight = "sunlight")
meeus_dn.df## # A tibble: 219,146 × 12
## day tz twilight.rise twilight.set longitude latitude
## <dttm> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1700-01-01 00:00:00 UTC -0.833 -0.833 0 51.5
## 2 1700-01-02 00:00:00 UTC -0.833 -0.833 0 51.5
## 3 1700-01-03 00:00:00 UTC -0.833 -0.833 0 51.5
## # ℹ 219,143 more rows
## # ℹ 6 more variables: address <chr>, sunrise <dttm>, noon <dttm>,
## # sunset <dttm>, daylength <dbl>, nightlength <dbl>Now we use the ‘suncalc’ package to compute the same values.
suncalc_dn.df <- getSunlightTimes(date = as.Date(times, tz = "UTC"), lat = 51.5, lon = 0)
head(suncalc_dn.df, n = 3L)## date lat lon solarNoon nadir
## 1 1700-01-01 51.5 0 1700-01-01 12:07:07 1700-01-01 00:07:07
## 2 1700-01-02 51.5 0 1700-01-02 12:07:32 1700-01-02 00:07:32
## 3 1700-01-03 51.5 0 1700-01-03 12:07:57 1700-01-03 00:07:57
## sunrise sunset sunriseEnd
## 1 1700-01-01 08:05:28 1700-01-01 16:08:45 1700-01-01 08:09:45
## 2 1700-01-02 08:05:02 1700-01-02 16:10:02 1700-01-02 08:09:19
## 3 1700-01-03 08:04:33 1700-01-03 16:11:21 1700-01-03 08:08:49
## sunsetStart dawn dusk
## 1 1700-01-01 16:04:28 1700-01-01 07:25:59 1700-01-01 16:48:14
## 2 1700-01-02 16:05:46 1700-01-02 07:25:39 1700-01-02 16:49:25
## 3 1700-01-03 16:07:06 1700-01-03 07:25:16 1700-01-03 16:50:38
## nauticalDawn nauticalDusk nightEnd
## 1 1700-01-01 06:43:24 1700-01-01 17:30:49 1700-01-01 06:02:59
## 2 1700-01-02 06:43:09 1700-01-02 17:31:55 1700-01-02 06:02:46
## 3 1700-01-03 06:42:51 1700-01-03 17:33:03 1700-01-03 06:02:31
## night goldenHourEnd goldenHour
## 1 1700-01-01 18:11:14 1700-01-01 09:05:03 1700-01-01 15:09:10
## 2 1700-01-02 18:12:18 1700-01-02 09:04:22 1700-01-02 15:10:42
## 3 1700-01-03 18:13:23 1700-01-03 09:03:38 1700-01-03 15:12:16We compute the differences in sun elevations, and their range. We plot the differences summarized as 2D density plots.
noon_differences <- suncalc_dn.df$solarNoon - meeus_dn.df$noon
range(differences)## [1] -3.672678 3.069026
ggplot(data.frame(times, noon_differences), aes(times, noon_differences)) +
stat_hdr() +
stat_panel_counts(label.x = "centre") +
labs(x = "Time", y = "Difference in noon time (seconds)")
sunrise_differences <- suncalc_dn.df$sunrise - meeus_dn.df$sunrise
range(sunrise_differences)## Time differences in secs
## [1] -653.0310 798.5153
ggplot(data.frame(times, sunrise_differences), aes(times, sunrise_differences)) +
stat_hdr() +
stat_panel_counts(label.x = "centre") +
labs(x = "Time", y = "Difference in sunrise time (seconds)")
So do these shifts affect the day length estimates?
suncalc_dn.df$daylength <- suncalc_dn.df$sunset - suncalc_dn.df$sunrise
daylength_differences <- suncalc_dn.df$daylength - meeus_dn.df$daylength
range(daylength_differences)## Time differences in hours
## [1] -0.3456996 0.3447686
ggplot(data.frame(times, daylength_differences), aes(times, daylength_differences)) +
stat_hdr() +
stat_panel_counts(label.x = "centre") +
labs(x = "Time", y = "Difference in daylength (h)")
The differences seem to be consistently increasing as we move away from the first third of the 21st century, 2016 or 2000 seem to be the least different depending on the value estimated. So the more accuracy needed, the shorter the period during which ‘suncalc’ provides useful values. In addition a large part of the innaccuracy seems to caused by increasing drift or bias. In “SunCalcMeus” the innaccuracy is smaller, cyclic through the year and constant from year to year for the 600 years considered.
It must be remembered that we here tested computed values based on different algorithms, not the observed apaprent times of sunrise and sunset. These vary, especially at high latitudes as a consequence of variation in atmospheric refraction that depends on current weather conditions at the time of observation.
Conclusion
This tests confirm that, in many but not all situations, the computation intensive algorithm from Meeus as implemented in package ‘SunCalcMeeus’ provides significantly better estimates of sunset and sunrise times, and of the position of the sun than package ‘suncalc’. This is achieved at the expense three times longer computation time, that is in most cases tolerable. The improvement in accuracy of using ‘SunCalcMeeus’ vs. ‘suncalc’ increases at times and dates farther away from the present.
Tests for dates further into the past would be needed. Obtaining reference data for a comparison requires more effort. Testing other geographic locations, specially higher latitudes is more important, and is in my to do list.