Функция 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"
Сейчас 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
2 ответа
Там нет ничего плохого в коде ядра, найденном в связанном посте, если вход when дается в UTC. Путаница заключалась в том, что ОП вошел не так Time Zone на сайте для Sys.time() из 2016-09-08 09:09:05 CDT:
Используя http://www.esrl.noaa.gov/gmd/grad/solcalc/ я получаю азимут 124 и высоту 38, что я считаю правильным.
Правильный Time Zone для входа на веб-сайт NOAA -5 за CDT ( см. этот сайт), который дает:
призвание 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
чтобы получить тот же результат. Обратите внимание, что нам нужно указать смещение от 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
за 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
Мне кажется, что долгота - это долгота относительно вашей позиции. Я не эксперт по этой теме, но в некотором смысле имеет смысл, что долгота была бы такой, поскольку она не сильно влияет на положение Солнца. Человек на данной широте в данное местное время увидит солнце в том же положении на небе, что и кто-то на другой долготе, но на той же широте и местном времени (игнорируя осложнения часовых поясов, представляющих собой сгруппированные области с дискретными границами и Землей). движется вокруг солнца).
Может быть, я не прочитал вопрос или функцию достаточно хорошо, но неожиданно, что долгота ведет себя так.
edit: чтение ответа @ aichao показывает, почему установка широты на ноль и установка местного времени работает приблизительно. Однако я не думаю, что это будет очень точно.