太阳位置的 R 函数给出了意想不到的结果

我想计算给定时间、纬度和经度的太阳的位置。我在这里找到了这个很棒的问题和答案：一天中给定时间、纬度和经度的太阳位置 https://stackoverflow.com/questions/8708048/position-of-the-sun-given-time-of-day-latitude-and-longitude/8823846#8823846。但是，当我评估该函数时，我得到了错误的结果。鉴于答案的质量，我几乎可以肯定我这边有问题，但我提出这个问题作为尝试解决问题的记录。

为了方便起见，下面重印了该函数的代码：

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)
  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)
  )
}

如果我运行：

sunPosition(when = Sys.time(),lat = 43, long = -89)

结果是：

  elevation azimuthJ azimuthC
1 -24.56604 55.26111 55.26111

Sys.time() 给出：

> Sys.time()
[1] "2016-09-08 09:09:05 CDT"

现在是上午九点，太阳已经升起。使用http://www.esrl.noaa.gov/gmd/grad/solcalc/ http://www.esrl.noaa.gov/gmd/grad/solcalc/我得到的方位角为 124，仰角为 38，我认为这是正确的。

我认为这可能是代码的问题，但我也从上面的答案中测试了 Josh 的原始 sunPosition 函数并得到了相同的结果。我的下一个想法是我的时间或时区有问题。

按照上述问题测试冬至，仍然给出了他们发现的相同结果并且看起来正确：

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

time <- as.POSIXct("2012-12-22 12:00:00")

sunPosition(when = time, lat = testPts$lat, long = testPts$long)

elevation azimuthJ azimuthC
1  72.43112 359.0787 359.0787
2  69.56493 180.7965 180.7965
3  63.56539 180.6247 180.6247
4  25.56642 180.3083 180.3083

当我进行相同的测试，但更改经度（-89）时，我在中午得到负海拔。

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

time <- as.POSIXct("2012-12-22 12:00:00 CDT")

sunPosition(when = time, lat = testPts$lat, long = testPts$long)

      elevation azimuthJ azimuthC
1  16.060136563 107.3420 107.3420
2   2.387033692 113.3522 113.3522
3   0.001378426 113.4671 113.4671
4 -14.190786786 108.8866 108.8866

链接帖子中找到的核心代码没有任何问题if输入when以 UTC 给出。令人困惑的是OP输入了错误Time Zone在网站中Sys.time() of 2016-09-08 09:09:05 CDT:

Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ http://www.esrl.noaa.gov/gmd/grad/solcalc/我得到的方位角为 124，仰角为 38，我认为这是正确的。

正确的Time Zone进入 NOAA 网站的输入是-5 for CDT (看这个网站 https://www.timeanddate.com/time/zones/cdt）， 这使：

Calling sunPosition将时间调整为 UTC 会得到类似的结果：

sunPosition(when = "2016-09-08 14:09:05", format="%Y-%m-%d %H:%M:%S",lat = 43, long = -89)
##  elevation azimuthJ azimuthC
##1  28.08683 110.4915 110.4915

现在，代码不会执行到 UTC 的转换。一种方法是将第一行替换为sunPosition:

if(is.character(when)) when <- strptime(when, format)

with

if(is.character(when)) 
  when <- strptime(when, format, tz="UTC")
else
  when <- as.POSIXlt(when, tz="UTC")

我们现在可以打电话sunPosition with:

sunPosition(when = "2016-09-08 09:09:05-0500", format="%Y-%m-%d %H:%M:%S%z",lat = 43, long = -89)
##  elevation azimuthJ azimuthC
##1  28.08683 110.4915 110.4915

得到相同的结果。请注意，我们NEED TO指定字符串文字中与 UTC 的偏移量format (%z) 打电话时sunPosition这边走。

随着这个改变sunPosition可以调用Sys.time()（我在东海岸）：

Sys.time()
##[1] "2016-09-08 12:42:08 EDT"
sunPosition(Sys.time(),lat = 43, long = -89)
##  elevation azimuthJ azimuthC
##1  49.24068 152.1195 152.1195

与 NOAA 网站相符

for Time Zone = -4 for EDT.