R函数给出的太阳位置结果与预期不符

Question

我希望能够根据时间、纬度和经度计算出太阳的位置。我在这里找到了一个非常好的问题和答案：Position of the sun given time of day, latitude and longitude 。但是，当我评估该函数时，得到的结果是不正确的。考虑到答案的质量，我几乎可以确定问题出在我这里，但我提出这个问题是为了记录尝试解决问题的过程。

以下是该函数的代码，方便查看：

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"

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

我以为可能是代码的问题，但我测试了乔希原来的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

- Tedward

---

我刚刚检查了所有24个时区的时间戳“2016-09-08 09:09:05 CDT”，没有一个给出“azimuthJ”为“55.26”或“azimuthC”为“55.26”。实际上，这两个值总是不同的。所以我认为你的代码可能有bug。 - Tim Biegeleisen

嗯，知道了。你按照我粘贴的代码运行了吗？结果是相同的坏结果吗？ - Tedward

我可以得到与你类似的结果，@Tedward: sunPosition(when = as.POSIXct("2016-09-08 09:09:05 CDT"),lat = 43, long = -89) 得出的结果为 -24.43463 55.45186 55.45186。 - ddunn801

Answer 1

如果输入的"when"是UTC时间，则链接文章中的核心代码没有问题。混淆之处在于OP在网站上为“2016-09-08 09:09:05 CDT”的Sys.time()输入了错误的时区：

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

输入到NOAA网站上的正确时区是CDT的-5（参见此网站），结果如下：

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

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)

使用

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

我们现在可以使用以下方式调用sunPosition函数：

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

为了获得相同的结果。请注意，在调用sunPosition时，我们需要在字符串文字和format（%z）中指定与UTC的偏移量。

更改后，可以使用Sys.time()（我在东海岸）调用sunPosition：

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网站匹配

对于时区为-4的EDT。

- aichao

---

谢谢你的回答，aichao！正如你所正确指出的那样，我的关键困惑是没有转换为UTC。 - Tedward

Answer 2

我认为问题出在经度上。如果我将经度设置为0，纬度设置为我的纬度，时间设置为我的时间，则可以得到正确的仰角和方位角数值。

> time <- Sys.time()
> time

[1] "2016-09-08 12:07:35 CDT"

> sunPosition(when = time, lat = 43, long = 0)

 elevation azimuthJ azimuthC
1  52.36687 184.1056 184.1056

在我看来，经度是相对于你的位置而言的。虽然我不是这个话题的专家，但从某种程度上讲，经度会以这种方式存在是有道理的，因为它对太阳位置没有太大影响。在给定纬度和当地时间下，一个人所见到的太阳位置与另一个经度不同但纬度和当地时间相同的人所见到的太阳位置是一样的（忽略时区被分成离散边界的区域和地球绕着太阳转动的复杂情况）。

也许我没有仔细阅读问题或函数，但经度这样运作是出乎意料的。

---

编辑： 阅读@aichao的答案可以解释为什么将纬度设置为零并使时间为本地时间可以近似地工作。但是，我认为这不会非常准确。

- Tedward