太阳位置的 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 的仰角，我认为这是正确的。

我想这可能是代码的问题，但我也从上面的答案中测试了 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

Answer 1

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

Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ I get an azimuth of 124 and elevation of 38, which I think is correct.

正确的Time Zone输入 NOAA 网站是 -5对于 CDT ( see this website )，它给出:



调用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)

和

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

得到相同的结果。请注意，我们 需要在字符串文字和 format 中指定与 UTC 的偏移量( %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 网站匹配



对于 Time Zone = -4对于 EDT .