R功能为阳光的位置提供意想不到的结果
我想计算太阳给定时间，纬度和经度的位置。 我在这里发现了这个伟大的问题和答案：太阳在一天中的时间，纬度和经度的位置。 但是，当我评估函数时，我得到不正确的结果。 鉴于答案的质量，我几乎可以肯定，我的问题出现了问题，但我把这个问题作为试图解决问题的记录。
为便于说明，以下是重印后功能的代码：
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/我得到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
如果在UTC中给出输入when ，链接帖子中发现的核心代码没有任何问题。 令人困惑的是，OP在Sys.time()的网站上输入了错误的Time Zone 2016-09-08 09:09:05 CDT ：
使用http://www.esrl.noaa.gov/gmd/grad/solcalc/我得到124的方位角和38的高度，我认为这是正确的。
输入NOAA网站的正确Time Zone为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这种方式调用sunPosition时，我们需要在字符串文本中指定UTC与UTC的偏移量， format （ %z ）。
随着这一变化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网站相匹配
对于EDT ， Time Zone = -4 。
我认为这个问题与经度有关。 如果我将经度设置为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的答案显示了为什么将纬度设置为零并使时间本地时间大致工作。 不过，我认为这不会很准确。
链接地址: http://www.djcxy.com/p/36223.html上一篇: R function for position of sun giving unexpected results
下一篇: embedding a tweet