太阳位置的 R 函数给出了意想不到的结果
R function for position of sun giving unexpected results
我想计算给定时间、纬度和经度的太阳位置。我在这里找到了这个很棒的问题和答案：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
在链接 post 中找到的核心代码没有任何问题 如果 输入 when 以 UTC 格式给出。令人困惑的是，OP 在 2016-09-08 09:09:05 CDT 的 Sys.time() 网站中输入了错误的 Time Zone：
Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ I get an azimuth of 124 and elevation of 38, which I think is correct.
输入 NOAA 网站的正确 Time Zone 是 -5 for CDT (see this website)，它给出：
将时间调整为 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
得到相同的结果。请注意，我们 需要 在字符串文字和 format (%z) 中以这种方式调用 sunPosition 时指定与 UTC 的偏移量。
通过此更改 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.
我认为问题出在经度上。如果我将经度设置为 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 的回答说明了为什么将纬度设置为零并使当地时间近似工作。但是，我认为这不会很准确。
我想计算给定时间、纬度和经度的太阳位置。我在这里找到了这个很棒的问题和答案：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
在链接 post 中找到的核心代码没有任何问题 如果 输入 when 以 UTC 格式给出。令人困惑的是，OP 在 2016-09-08 09:09:05 CDT 的 Sys.time() 网站中输入了错误的 Time Zone：
Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ I get an azimuth of 124 and elevation of 38, which I think is correct.
输入 NOAA 网站的正确 Time Zone 是 -5 for CDT (see this website)，它给出：
将时间调整为 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
得到相同的结果。请注意，我们 需要 在字符串文字和 format (%z) 中以这种方式调用 sunPosition 时指定与 UTC 的偏移量。
通过此更改 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.
我认为问题出在经度上。如果我将经度设置为 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 的回答说明了为什么将纬度设置为零并使当地时间近似工作。但是，我认为这不会很准确。