R功能为太阳位置带来意想不到的效果

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.

现在是上午9点,太阳升了。使用http://www.esrl.noaa.gov/gmd/grad/solcalc/我的方位角为124,高程为38,我认为这是正确的。

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.

我想也许这是代码的问题,但我也从上面的答案中测试了Josh的原始sunPosition函数并得到了相同的结果。我的下一个想法是我的时间或时区存在问题。

---

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.

当我做同样的测试,但改变经度(-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 个解决方案

#1

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:

如果以UTC格式给出输入,则链接帖子中找到的核心代码没有任何问题。令人困惑的是,OP在2016-09-08 09:09:05 CDT的Sys.time()网站进入了错误的时区:

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

使用http://www.esrl.noaa.gov/gmd/grad/solcalc/我的方位角为124,高程为38,我认为这是正确的。

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

输入NOAA网站的正确时区为-5(CDT)(见本网站),其中给出:

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

在调整为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

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

现在,代码不会将此转换为UTC。一种方法是替换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:

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.

得到相同的结果。请注意,我们需要在以这种方式调用sunPosition时以字符串文字和格式(%z)指定UTC的偏移量。

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

通过此更改,可以使用Sys.time()(我在东海岸)调用sunPosition:

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

与NOAA网站相匹配

for Time Zone = -4 for EDT.

对于EDT,时区= -4。

#2

I think the issue is with longitude. If I set longitude to 0 and latitude to my latitude and time to my time, I get correct values for elevation and azimuth.

我认为这个问题是经度问题。如果我将经度设置为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

It seems to me that longitude is longitude relative to your position. I'm not an expert on this topic, but in a way it makes sense that longitude would be this way since it doesn't have much effect on solar position. A person at a given latitude at a given local time will see the sun in the same position in the sky as someone else at a different longitude, but same latitude and local time (ignoring complications of timezones being binned regions with discrete boundaries and the earth moving around the sun).

在我看来,经度是相对于你的位置的经度。我不是这个主题的专家,但从某种意义上说,经度是这样的,因为它对太阳位置没有太大影响。在给定的当地时间给定纬度的人将看到太阳与天空中的其他人在不同的经度上相同的位置,但相同的纬度和当地时间(忽略时区的复杂性是具有离散边界和地球的分区区域在太阳周围移动)。

Maybe I didn't read the question or the function well enough, but it is unexpected that longitude behaves this way.

也许我没有很好地阅读这个问题或功能,但经度表现出来是出乎意料的。

---

edit: Reading @aichao's answer shows why setting latitude to zero and making the time local time works approximately. However, I don't think this will be very accurate.

编辑:阅读@ aichao的答案显示为什么将纬度设置为零并使当地时间大致工作。但是,我认为这不会很准确。