小开

最佳答案

This seems like an important topic, so I've posted a longer than typical answer: if this algorithm is to be used by others in the future, I think it's important that it be accompanied by references to the literature from which it has been derived.

The short answer

As you've noted, your posted code does not work properly for locations near the equator, or in the southern hemisphere.

To fix it, simply replace these lines in your original code:

elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

with these:

cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]

It should now work for any location on the globe.

Discussion

The code in your example is adapted almost verbatim from a 1988 article by J.J. Michalsky (Solar Energy. 40:227-235). That article in turn refined an algorithm presented in a 1978 article by R. Walraven (Solar Energy. 20:393-397). Walraven reported that the method had been used successfully for several years to precisely position a polarizing radiometer in Davis, CA (38° 33' 14" N, 121° 44' 17" W).

Both Michalsky's and Walraven's code contains important/fatal errors. In particular, while Michalsky's algorithm works just fine in most of the United States, it fails (as you've found) for areas near the equator, or in the southern hemisphere. In 1989, J.W. Spencer of Victoria, Australia, noted the same thing (Solar Energy. 42(4):353):

Dear Sir:

Michalsky's method for assigning the calculated azimuth to the correct quadrant, derived from Walraven, does not give correct values when applied for Southern (negative) latitudes. Further the calculation of the critical elevation (elc) will fail for a latitude of zero because of division by zero. Both these objections can be avoided simply by assigning the azimuth to the correct quadrant by considering the sign of cos(azimuth).

My edits to your code are based on the corrections suggested by Spencer in that published Comment. I have simply altered them somewhat to ensure that the R function sunPosition() remains 'vectorized' (i.e. working properly on vectors of point locations, rather than needing to be passed one point at a time).

Accuracy of the function `sunPosition()`

To test that sunPosition() works correctly, I've compared its results with those calculated by the National Oceanic and Atmospheric Administration's Solar Calculator. In both cases, sun positions were calculated for midday (12:00 PM) on the southern summer solstice (December 22nd), 2012. All results were in agreement to within 0.02 degrees.

testPts <- data.frame(lat = c(-41,-3,3, 41),
long = c(0, 0, 0, 0))


# Sun's position as returned by the NOAA Solar Calculator,
NOAA <- data.frame(elevNOAA = c(72.44, 69.57, 63.57, 25.6),
azNOAA = c(359.09, 180.79, 180.62, 180.3))


# Sun's position as returned by sunPosition()
sunPos <- sunPosition(year = 2012,
month = 12,
day = 22,
hour = 12,
min = 0,
sec = 0,
lat = testPts$lat,
long = testPts$long)


cbind(testPts, NOAA, sunPos)
#   lat long elevNOAA azNOAA elevation  azimuth
# 1 -41    0    72.44 359.09  72.43112 359.0787
# 2  -3    0    69.57 180.79  69.56493 180.7965
# 3   3    0    63.57 180.62  63.56539 180.6247
# 4  41    0    25.60 180.30  25.56642 180.3083

Other errors in the code

There are at least two other (quite minor) errors in the posted code. The first causes February 29th and March 1st of leap years to both be tallied as day 61 of the year. The second error derives from a typo in the original article, which was corrected by Michalsky in a 1989 note (Solar Energy. 43(5):323).

This code block shows the offending lines, commented out and followed immediately by corrected versions:

# leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
day >= 60 & !(month==2 & day==60)


# oblqec <- 23.429 - 0.0000004 * time
oblqec <- 23.439 - 0.0000004 * time

Corrected version of `sunPosition()`

Here is the corrected code that was verified above:

sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
lat=46.5, long=6.5) {


twopi <- 2 * pi
deg2rad <- pi / 180


# Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
day <- day + cumsum(month.days)[month]
leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) &
day >= 60 & !(month==2 & day==60)
day[leapdays] <- day[leapdays] + 1


# Get Julian date - 2400000
hour <- hour + min / 60 + sec / 3600 # hour plus fraction
delta <- year - 1949
leap <- trunc(delta / 4) # former leapyears
jd <- 32916.5 + delta * 365 + leap + day + hour / 24


# The input to the Atronomer's almanach is the difference between
# the Julian date and JD 2451545.0 (noon, 1 January 2000)
time <- jd - 51545.


# Ecliptic coordinates


# Mean longitude
mnlong <- 280.460 + .9856474 * time
mnlong <- mnlong %% 360
mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360


# Mean anomaly
mnanom <- 357.528 + .9856003 * time
mnanom <- mnanom %% 360
mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
mnanom <- mnanom * deg2rad


# Ecliptic longitude and obliquity of ecliptic
eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
eclong <- eclong %% 360
eclong[eclong < 0] <- eclong[eclong < 0] + 360
oblqec <- 23.439 - 0.0000004 * time
eclong <- eclong * deg2rad
oblqec <- oblqec * deg2rad


# Celestial coordinates
# Right ascension and declination
num <- cos(oblqec) * sin(eclong)
den <- cos(eclong)
ra <- atan(num / den)
ra[den < 0] <- ra[den < 0] + pi
ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + twopi
dec <- asin(sin(oblqec) * sin(eclong))


# Local coordinates
# Greenwich mean sidereal time
gmst <- 6.697375 + .0657098242 * time + hour
gmst <- gmst %% 24
gmst[gmst < 0] <- gmst[gmst < 0] + 24.


# Local mean sidereal time
lmst <- gmst + long / 15.
lmst <- lmst %% 24.
lmst[lmst < 0] <- lmst[lmst < 0] + 24.
lmst <- lmst * 15. * deg2rad


# Hour angle
ha <- lmst - ra
ha[ha < -pi] <- ha[ha < -pi] + twopi
ha[ha > pi] <- ha[ha > pi] - twopi


# Latitude to radians
lat <- lat * deg2rad


# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
az <- asin(-cos(dec) * sin(ha) / cos(el))


# For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]


# if (0 < sin(dec) - sin(el) * sin(lat)) {
#     if(sin(az) < 0) az <- az + twopi
# } else {
#     az <- pi - az
# }




el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad


return(list(elevation=el, azimuth=az))
}

References:

Michalsky, J.J. 1988. The Astronomical Almanac's algorithm for approximate solar position (1950-2050). Solar Energy. 40(3):227-235.

Michalsky, J.J. 1989. Errata. Solar Energy. 43(5):323.

Spencer, J.W. 1989. Comments on "The Astronomical Almanac's Algorithm for Approximate Solar Position (1950-2050)". Solar Energy. 42(4):353.

Walraven, R. 1978. Calculating the position of the sun. Solar Energy. 20:393-397.

小开

Using "NOAA Solar Calculations" from one of the links above I have changed a bit the final part of the function by using a slighly different algorithm that, I hope, have translated without errors. I have commented out the now-useless code and added the new algorithm just after the latitude to radians conversion:

# -----------------------------------------------
# New code
# Solar zenith angle
zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
# Solar azimuth
az <- acos(((sin(lat) * cos(zenithAngle)) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
rm(zenithAngle)
# -----------------------------------------------


# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
#az <- asin(-cos(dec) * sin(ha) / cos(el))
#elc <- asin(sin(dec) / sin(lat))
#az[el >= elc] <- pi - az[el >= elc]
#az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi


el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad


# -----------------------------------------------
# New code
if (ha > 0) az <- az + 180 else az <- 540 - az
az <- az %% 360
# -----------------------------------------------


return(list(elevation=el, azimuth=az))

To verify azimuth trend in the four cases you mentioned let's plot it against time of day:

hour <- seq(from = 0, to = 23, by = 0.5)
azimuth <- data.frame(hour = hour)
az41S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-41,0)$azimuth)
az03S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-03,0)$azimuth)
az03N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,03,0)$azimuth)
az41N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,41,0)$azimuth)
azimuth <- cbind(azimuth, az41S, az03S, az41N, az03N)
rm(az41S, az03S, az41N, az03N)
library(ggplot2)
azimuth.plot <- melt(data = azimuth, id.vars = "hour")
ggplot(aes(x = hour, y = value, color = variable), data = azimuth.plot) +
geom_line(size = 2) +
geom_vline(xintercept = 12) +
facet_wrap(~ variable)

Image attached:

enter image description here

小开

This is a suggested update to Josh's excellent answer.

Much of the start of the function is boilerplate code for calculating the number of days since midday on 1st Jan 2000. This is much better dealt with using R's existing date and time function.

I also think that rather than having six different variables to specify the date and time, it's easier (and more consistent with other R functions) to specify an existing date object or a date strings + format strings.

Here are two helper functions

astronomers_almanac_time <- function(x)
{
origin <- as.POSIXct("2000-01-01 12:00:00")
as.numeric(difftime(x, origin, units = "days"))
}


hour_of_day <- function(x)
{
x <- as.POSIXlt(x)
with(x, hour + min / 60 + sec / 3600)
}

And the start of the function now simplifies to

sunPosition <- function(when = Sys.time(), format, lat=46.5, long=6.5) {


twopi <- 2 * pi
deg2rad <- pi / 180


if(is.character(when)) when <- strptime(when, format)
time <- astronomers_almanac_time(when)
hour <- hour_of_day(when)
#...

The other oddity is in the lines like

mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

Since mnlong has had %% called on its values, they should all be non-negative already, so this line is superfluous.

小开

Here's a rewrite in that's more idiomatic to R, and easier to debug and maintain. It is essentially Josh's answer, but with azimuth calculated using both Josh and Charlie's algorithms for comparison. I've also included the simplifications to the date code from my other answer. The basic principle was to split the code up into lots of smaller functions that you can more easily write unit tests for.

astronomersAlmanacTime <- function(x)
{
# Astronomer's almanach time is the number of
# days since (noon, 1 January 2000)
origin <- as.POSIXct("2000-01-01 12:00:00")
as.numeric(difftime(x, origin, units = "days"))
}


hourOfDay <- function(x)
{
x <- as.POSIXlt(x)
with(x, hour + min / 60 + sec / 3600)
}


degreesToRadians <- function(degrees)
{
degrees * pi / 180
}


radiansToDegrees <- function(radians)
{
radians * 180 / pi
}


meanLongitudeDegrees <- function(time)
{
(280.460 + 0.9856474 * time) %% 360
}


meanAnomalyRadians <- function(time)
{
degreesToRadians((357.528 + 0.9856003 * time) %% 360)
}


eclipticLongitudeRadians <- function(mnlong, mnanom)
{
degreesToRadians(
(mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)) %% 360
)
}


eclipticObliquityRadians <- function(time)
{
degreesToRadians(23.439 - 0.0000004 * time)
}


rightAscensionRadians <- function(oblqec, eclong)
{
num <- cos(oblqec) * sin(eclong)
den <- cos(eclong)
ra <- atan(num / den)
ra[den < 0] <- ra[den < 0] + pi
ra[den >= 0 & num < 0] <- ra[den >= 0 & num < 0] + 2 * pi
ra
}


rightDeclinationRadians <- function(oblqec, eclong)
{
asin(sin(oblqec) * sin(eclong))
}


greenwichMeanSiderealTimeHours <- function(time, hour)
{
(6.697375 + 0.0657098242 * time + hour) %% 24
}


localMeanSiderealTimeRadians <- function(gmst, long)
{
degreesToRadians(15 * ((gmst + long / 15) %% 24))
}


hourAngleRadians <- function(lmst, ra)
{
((lmst - ra + pi) %% (2 * pi)) - pi
}


elevationRadians <- function(lat, dec, ha)
{
asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
}


solarAzimuthRadiansJosh <- function(lat, dec, ha, el)
{
az <- asin(-cos(dec) * sin(ha) / cos(el))
cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + 2 * pi
az[!cosAzPos] <- pi - az[!cosAzPos]
az
}


solarAzimuthRadiansCharlie <- function(lat, dec, ha)
{
zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
az <- acos((sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
ifelse(ha > 0, az + pi, 3 * pi - az) %% (2 * pi)
}


sunPosition <- function(when = Sys.time(), format, lat = 46.5, long = 6.5)
{
if(is.character(when)) when <- strptime(when, format)
when <- lubridate::with_tz(when, "UTC")
time <- astronomersAlmanacTime(when)
hour <- hourOfDay(when)


# Ecliptic coordinates
mnlong <- meanLongitudeDegrees(time)
mnanom <- meanAnomalyRadians(time)
eclong <- eclipticLongitudeRadians(mnlong, mnanom)
oblqec <- eclipticObliquityRadians(time)


# Celestial coordinates
ra <- rightAscensionRadians(oblqec, eclong)
dec <- rightDeclinationRadians(oblqec, eclong)


# Local coordinates
gmst <- greenwichMeanSiderealTimeHours(time, hour)
lmst <- localMeanSiderealTimeRadians(gmst, long)


# Hour angle
ha <- hourAngleRadians(lmst, ra)


# Latitude to radians
lat <- degreesToRadians(lat)


# Azimuth and elevation
el <- elevationRadians(lat, dec, ha)
azJ <- solarAzimuthRadiansJosh(lat, dec, ha, el)
azC <- solarAzimuthRadiansCharlie(lat, dec, ha)


data.frame(
elevation = radiansToDegrees(el),
azimuthJ  = radiansToDegrees(azJ),
azimuthC  = radiansToDegrees(azC)
)
}

小开

I needed sun position in a Python project. I adapted Josh O'Brien's algorithm.

Thank you Josh.

In case it could be useful to anyone, here's my adaptation.

Note that my project only needed instant sun position so time is not a parameter.

def sunPosition(lat=46.5, long=6.5):


# Latitude [rad]
lat_rad = math.radians(lat)


# Get Julian date - 2400000
day = time.gmtime().tm_yday
hour = time.gmtime().tm_hour + \
time.gmtime().tm_min/60.0 + \
time.gmtime().tm_sec/3600.0
delta = time.gmtime().tm_year - 1949
leap = delta / 4
jd = 32916.5 + delta * 365 + leap + day + hour / 24


# The input to the Atronomer's almanach is the difference between
# the Julian date and JD 2451545.0 (noon, 1 January 2000)
t = jd - 51545


# Ecliptic coordinates


# Mean longitude
mnlong_deg = (280.460 + .9856474 * t) % 360


# Mean anomaly
mnanom_rad = math.radians((357.528 + .9856003 * t) % 360)


# Ecliptic longitude and obliquity of ecliptic
eclong = math.radians((mnlong_deg +
1.915 * math.sin(mnanom_rad) +
0.020 * math.sin(2 * mnanom_rad)
) % 360)
oblqec_rad = math.radians(23.439 - 0.0000004 * t)


# Celestial coordinates
# Right ascension and declination
num = math.cos(oblqec_rad) * math.sin(eclong)
den = math.cos(eclong)
ra_rad = math.atan(num / den)
if den < 0:
ra_rad = ra_rad + math.pi
elif num < 0:
ra_rad = ra_rad + 2 * math.pi
dec_rad = math.asin(math.sin(oblqec_rad) * math.sin(eclong))


# Local coordinates
# Greenwich mean sidereal time
gmst = (6.697375 + .0657098242 * t + hour) % 24
# Local mean sidereal time
lmst = (gmst + long / 15) % 24
lmst_rad = math.radians(15 * lmst)


# Hour angle (rad)
ha_rad = (lmst_rad - ra_rad) % (2 * math.pi)


# Elevation
el_rad = math.asin(
math.sin(dec_rad) * math.sin(lat_rad) + \
math.cos(dec_rad) * math.cos(lat_rad) * math.cos(ha_rad))


# Azimuth
az_rad = math.asin(
- math.cos(dec_rad) * math.sin(ha_rad) / math.cos(el_rad))


if (math.sin(dec_rad) - math.sin(el_rad) * math.sin(lat_rad) < 0):
az_rad = math.pi - az_rad
elif (math.sin(az_rad) < 0):
az_rad += 2 * math.pi


return el_rad, az_rad

小开

I encountered a slight problem with a data point & Richie Cotton's functions above (in the implementation of Charlie's code)

longitude= 176.0433687000000020361767383292317390441894531250
latitude= -39.173830619999996827118593500927090644836425781250
event_time = as.POSIXct("2013-10-24 12:00:00", format="%Y-%m-%d %H:%M:%S", tz = "UTC")
sunPosition(when=event_time, lat = latitude, long = longitude)
elevation azimuthJ azimuthC
1 -38.92275      180      NaN
Warning message:
In acos((sin(lat) * cos(zenithAngle) - sin(dec))/(cos(lat) * sin(zenithAngle))) : NaNs produced

because in the solarAzimuthRadiansCharlie function there has been floating point excitement around an angle of 180 such that (sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle)) is the tiniest amount over 1, 1.0000000000000004440892098, which generates a NaN as the input to acos should not be above 1 or below -1.

I suspect there might be similar edge cases for Josh's calculation, where floating point rounding effects cause the input for the asin step to be outside -1:1 but I have not hit them in my particular dataset.

In the half-dozen or so cases I have hit this, the "true" (middle of the day or night) is when the issue occurs so empirically the true value should be 1/-1. For that reason, I would be comfortable fixing that by applying a rounding step within solarAzimuthRadiansJosh and solarAzimuthRadiansCharlie. I'm not sure what the theoretical accuracy of the NOAA algorithm is (the point at which numerical accuracy stops mattering anyway) but rounding to 12 decimal places fixed the data in my data set.

给定时间的太阳位置，经纬度