R功能为阳光的位置提供意想不到的结果 [英] R function for position of sun giving unexpected results

问题描述

我想计算太阳给定时间，纬度和经度的位置。我在这里找到了这个很好的问题和答案：），它给出：

调用时间调整为UTC时调用 sunPosition 给出了类似的结果： / p>

  sunPosition（when =2016-09-08 14:09:05，format =％Y-％m-％ d％H：％M：％S，lat = 43，long = -89）
 ##仰角azimuthJ方位角C 
 ## 1 28.08683 110.4915 110.4915 
  

现在，代码不会执行UTC转换。一种方法是替换 sunPosition 中的第一行：

  if（is.character（when））when<  -  strptime（when，format）
  

with

$ b $

  if（is.character（when））
当<  -  strptime（when，format，tz = UTC）
 else 
当<  -  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）
 ##提升azimuthJ azimuthC 
 ## 1 28.08683 110.4915 110.4915 
  

得到相同的结果。请注意，我们 NEED TO 指定字符串文本和格式中UTC的偏移量（％z sunPosition 以这种方式调用时，c $ c>）。



通过这个改变 sunPosition 可以用 Sys.time（）（我在东海岸）调用：



  Sys.time（）
 ## [1]2016-09-08 12:42:08 EDT
 sunPosition（Sys.time（） ，lat = 43，long = -89）
 ##仰角azimuthJ azimuthC 
 ## 1 49.24068 152.1195 152.1195 
  

与NOAA网站匹配



代表时区 = -4 代表 EDT 。

I'd like to calculate the position of the sun given time, latitude, and longitude. I found this great question and answer here: Position of the sun given time of day, latitude and longitude. However, when I evaluate the function I get incorrect results. Given the quality of the answer, I'm almost certain there's something wrong on my end, but I'm asking this question as a record of trying to solve the problem.

Here's the code for the function reprinted below for convenience:

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)
  )
}

If I run:

sunPosition(when = Sys.time(),lat = 43, long = -89)

The result is:

  elevation azimuthJ azimuthC
1 -24.56604 55.26111 55.26111

Sys.time() gives:

> Sys.time()
[1] "2016-09-08 09:09:05 CDT"

It's 9am, and the sun is up. Using http://www.esrl.noaa.gov/gmd/grad/solcalc/ I get an azimuth of 124 and elevation of 38, which I think is correct.

I thought maybe it was an issue with the code, but I also tested Josh's original sunPosition function from the above answer and got the same results. My next thought is that there is an issue with my time or timezone.

---

testing the winter solstice as done in the above question, still gives the same results they found and looks correct:

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

When I do the same test, but change the longitude (-89), I get a negative elevation at noon.

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

解决方案

There is nothing wrong with the core code found in the linked post if the input when is given in UTC. The confusion was that the OP entered the wrong Time Zone in the website for the Sys.time() of 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.

The correct Time Zone to input into the NOAA website is -5 for CDT (see this website), which gives:

Calling sunPosition with the time adjusted to UTC gives a similar result:

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

Now, the code does not do this conversion to UTC. One way to do that is to replace the first line in sunPosition:

if(is.character(when)) when <- strptime(when, format)

with

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

We can now call sunPosition with:

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

to get the same result. Note that we NEED TO specify the offset from UTC in the string literal and in the format (%z) when calling sunPosition this way.

With this change sunPosition can be called with Sys.time() (I'm on the East Coast):

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

which matches the NOAA website

for Time Zone = -4 for EDT.

这篇关于R功能为阳光的位置提供意想不到的结果的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！