r - 經線英文 - 太陽的位置,一天中的時間,經度和緯度



經緯度0 0 (4)

這個問題在三年前已經問before了。 有一個答案,但我發現在解決方案中有一個小故障。

下面的代碼在R.我已經將它移植到另一種語言,但是已經直接在R中測試了原始代碼,以確保問題不在我的移植中。

sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
                    lat=46.5, long=6.5) {


  twopi <- 2 * pi
  deg2rad <- pi / 180

  # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
  month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
  day <- day + cumsum(month.days)[month]
  leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
  day[leapdays] <- day[leapdays] + 1

  # Get Julian date - 2400000
  hour <- hour + min / 60 + sec / 3600 # hour plus fraction
  delta <- year - 1949
  leap <- trunc(delta / 4) # former leapyears
  jd <- 32916.5 + delta * 365 + leap + day + hour / 24

  # The input to the Atronomer's almanach is the difference between
  # the Julian date and JD 2451545.0 (noon, 1 January 2000)
  time <- jd - 51545.

  # Ecliptic coordinates

  # Mean longitude
  mnlong <- 280.460 + .9856474 * time
  mnlong <- mnlong %% 360
  mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

  # Mean anomaly
  mnanom <- 357.528 + .9856003 * time
  mnanom <- mnanom %% 360
  mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
  mnanom <- mnanom * deg2rad

  # Ecliptic longitude and obliquity of ecliptic
  eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
  eclong <- eclong %% 360
  eclong[eclong < 0] <- eclong[eclong < 0] + 360
  oblqec <- 23.429 - 0.0000004 * time
  eclong <- eclong * deg2rad
  oblqec <- oblqec * deg2rad

  # Celestial coordinates
  # Right ascension and declination
  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] + twopi
  dec <- asin(sin(oblqec) * sin(eclong))

  # Local coordinates
  # Greenwich mean sidereal time
  gmst <- 6.697375 + .0657098242 * time + hour
  gmst <- gmst %% 24
  gmst[gmst < 0] <- gmst[gmst < 0] + 24.

  # Local mean sidereal time
  lmst <- gmst + long / 15.
  lmst <- lmst %% 24.
  lmst[lmst < 0] <- lmst[lmst < 0] + 24.
  lmst <- lmst * 15. * deg2rad

  # Hour angle
  ha <- lmst - ra
  ha[ha < -pi] <- ha[ha < -pi] + twopi
  ha[ha > pi] <- ha[ha > pi] - twopi

  # Latitude to radians
  lat <- lat * deg2rad

  # Azimuth and elevation
  el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
  az <- asin(-cos(dec) * sin(ha) / cos(el))
  elc <- asin(sin(dec) / sin(lat))
  az[el >= elc] <- pi - az[el >= elc]
  az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

  el <- el / deg2rad
  az <- az / deg2rad
  lat <- lat / deg2rad

  return(list(elevation=el, azimuth=az))
}

我遇到的問題是它返回的方位角似乎是錯誤的。 例如,如果我在(南部)夏至12點在0ºE和41ºS,3ºS,3ºN和41ºN位置運行函數:

> sunPosition(2012,12,22,12,0,0,-41,0)
$elevation
[1] 72.42113

$azimuth
[1] 180.9211

> sunPosition(2012,12,22,12,0,0,-3,0)
$elevation
[1] 69.57493

$azimuth
[1] -0.79713

Warning message:
In asin(sin(dec)/sin(lat)) : NaNs produced
> sunPosition(2012,12,22,12,0,0,3,0)
$elevation
[1] 63.57538

$azimuth
[1] -0.6250971

Warning message:
In asin(sin(dec)/sin(lat)) : NaNs produced
> sunPosition(2012,12,22,12,0,0,41,0)
$elevation
[1] 25.57642

$azimuth
[1] 180.3084

這些數字似乎並不正確。 我很滿意的海拔 - 前兩個應該大致相同,第三個降低,第四個降低。 然而,第一個方位角應該大致為北,而它給出的數字是完全相反的。 剩下的三個應該大致指向南方,但只有最後一個。 中間的兩個點位於北部,再次是180度。

正如你所看到的,還有一些低緯度地區引發的誤差(靠近赤道)

我相信這部分錯誤在第三行(以elc )觸發。

  # Azimuth and elevation
  el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
  az <- asin(-cos(dec) * sin(ha) / cos(el))
  elc <- asin(sin(dec) / sin(lat))
  az[el >= elc] <- pi - az[el >= elc]
  az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

我搜索了一下,並在C中找到了類似的代碼塊,將其轉換為R它用來計算方位角的線將是類似於

az <- atan(sin(ha) / (cos(ha) * sin(lat) - tan(dec) * cos(lat)))

這裡的輸出似乎正朝著正確的方向前進,但當我們將其轉化為度數時,我無法一直給予正確的答案。

對代碼進行修正(懷疑它只是上面的幾行),使其計算出正確的方位角會非常棒。

https://src-bin.com


Answer #1

使用上述鏈接中的“NOAA太陽能計算”,我已經通過使用一種略微不同的算法改變了函數的最後部分,我希望這些算法能夠毫無差錯地進行翻譯。 我已經將現在無用的代碼註釋掉了,並在緯度經過弧度轉換之後添加了新算法:

# -----------------------------------------------
# New code
# Solar zenith angle
zenithAngle <- acos(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
# Solar azimuth
az <- acos(((sin(lat) * cos(zenithAngle)) - sin(dec)) / (cos(lat) * sin(zenithAngle)))
rm(zenithAngle)
# -----------------------------------------------

# Azimuth and elevation
el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
#az <- asin(-cos(dec) * sin(ha) / cos(el))
#elc <- asin(sin(dec) / sin(lat))
#az[el >= elc] <- pi - az[el >= elc]
#az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

el <- el / deg2rad
az <- az / deg2rad
lat <- lat / deg2rad

# -----------------------------------------------
# New code
if (ha > 0) az <- az + 180 else az <- 540 - az
az <- az %% 360
# -----------------------------------------------

return(list(elevation=el, azimuth=az))

為了驗證您提到的四種情況下的方位角趨勢,我們將其與一天中的時間進行對比:

hour <- seq(from = 0, to = 23, by = 0.5)
azimuth <- data.frame(hour = hour)
az41S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-41,0)$azimuth)
az03S <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,-03,0)$azimuth)
az03N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,03,0)$azimuth)
az41N <- apply(azimuth, 1, function(x) sunPosition(2012,12,22,x,0,0,41,0)$azimuth)
azimuth <- cbind(azimuth, az41S, az03S, az41N, az03N)
rm(az41S, az03S, az41N, az03N)
library(ggplot2)
azimuth.plot <- melt(data = azimuth, id.vars = "hour")
ggplot(aes(x = hour, y = value, color = variable), data = azimuth.plot) + 
    geom_line(size = 2) + 
    geom_vline(xintercept = 12) + 
    facet_wrap(~ variable)

附圖:


Answer #2

我在數據點和Richie Cotton的函數上面遇到了一個小問題(在執行Charlie的代碼時)

longitude= 176.0433687000000020361767383292317390441894531250
latitude= -39.173830619999996827118593500927090644836425781250
event_time = as.POSIXct("2013-10-24 12:00:00", format="%Y-%m-%d %H:%M:%S", tz = "UTC")
sunPosition(when=event_time, lat = latitude, long = longitude)
elevation azimuthJ azimuthC
1 -38.92275      180      NaN
Warning message:
In acos((sin(lat) * cos(zenithAngle) - sin(dec))/(cos(lat) * sin(zenithAngle))) : NaNs produced

因為在太陽和ZimuthRadiansCharlie函數中,浮點激勵的角度約為180°,因此(sin(lat) * cos(zenithAngle) - sin(dec)) / (cos(lat) * sin(zenithAngle))是最小的量超過1,1.0000000000000004440892098,它產生一個NaN作為acos的輸入不應該高於1或低於-1。

我懷疑Josh的計算可能會有類似的邊緣情況,其中浮點舍入效應會導致asin步驟的輸入超出-1:1,但我沒有在我的特定數據集中擊中它們。

在我碰到這種情況時,“真實的”(白天或黑夜中)是這樣的情況,當問題發生如此經驗時,真實值應該是1 / -1。 出於這個原因,我很樂意通過在solarAzimuthRadiansJoshsolarAzimuthRadiansCharlie採用舍入步驟來solarAzimuthRadiansJosh solarAzimuthRadiansCharlie 。 我不確定NOAA算法的理論精度是多少(數值精度無論如何停止的點),但舍入到小數點後12位將數據固定在我的數據集中。


Answer #3

這似乎是一個重要的話題,所以我發布了一個比典型的答案更長的時間:如果這個算法將來會被其他人使用,我認為重要的是它伴隨著對其衍生出的文獻的參考。

簡短的答案

如您所知,您的發布代碼在赤道附近或南半球的位置無法正常工作。

要修復它,只需在原始代碼中替換這些行即可:

elc <- asin(sin(dec) / sin(lat))
az[el >= elc] <- pi - az[el >= elc]
az[el <= elc & ha > 0] <- az[el <= elc & ha > 0] + twopi

用這些:

cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
sinAzNeg <- (sin(az) < 0)
az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
az[!cosAzPos] <- pi - az[!cosAzPos]

它現在應該適用於全球任何地點。

討論

在你的例子中的代碼幾乎逐字從1988年JJ Michalsky的文章(太陽能40:227-235)改編而來。 該文章反過來改進了R. Walraven在1978年的一篇文章中提出的算法(Solar Energy 20:393-397)。 Walraven報導說,這種方法已經成功使用了好幾年,可以在加利福尼亞州戴維斯(38°33'14“N,121°44'17”W)精確定位偏振輻射計。

Michalsky和Walraven的代碼都包含重要的/致命的錯誤。 特別是,雖然Michalsky的算法在美國大部分地區都可以正常工作,但它在赤道附近或南半球的地區失敗了(如您所見)。 1989年,澳大利亞維多利亞的JW Spencer指出了同樣的事情(太陽能42(4):353):

親愛的先生:

Michalsky用於將計算出的方位角分配給來自Walraven的正確象限的方法在應用於南方(負)緯度時沒有給出正確的值。 此外,臨界高程(elc)的計算由於被零除而在緯度為零時將失敗。 通過考慮cos(方位角)的符號,通過將方位角分配給正確的象限,可以避免這兩種異議。

我對您的代碼的編輯基於Spencer在該已發布的評論中提出的更正。 為了確保R函數sunPosition()保持“向量化”(即在點位置的向量上正常工作,而不是一次只傳遞一個點),我只是稍微改變了它們。

函數sunPosition()精度

為了測試sunPosition()是否正常工作,我將其結果與國家海洋和大氣管理局的太陽能計算器計算結果進行了比較 。 在這兩種情況下,2012年南部夏至(12月22日)的中午(中午12點)的日照位置都計算在內。所有結果都一致在0.02度以內。

testPts <- data.frame(lat = c(-41,-3,3, 41), 
                      long = c(0, 0, 0, 0))

# Sun's position as returned by the NOAA Solar Calculator,
NOAA <- data.frame(elevNOAA = c(72.44, 69.57, 63.57, 25.6),
                   azNOAA = c(359.09, 180.79, 180.62, 180.3))

# Sun's position as returned by sunPosition()
sunPos <- sunPosition(year = 2012,
                      month = 12,
                      day = 22,
                      hour = 12,
                      min = 0,
                      sec = 0,
                      lat = testPts$lat,
                      long = testPts$long)

cbind(testPts, NOAA, sunPos)
#   lat long elevNOAA azNOAA elevation  azimuth
# 1 -41    0    72.44 359.09  72.43112 359.0787
# 2  -3    0    69.57 180.79  69.56493 180.7965
# 3   3    0    63.57 180.62  63.56539 180.6247
# 4  41    0    25.60 180.30  25.56642 180.3083

代碼中的其他錯誤

在發布的代碼中至少有兩個其他(很小的)錯誤。 第一個原因是2月29日和3月1日的閏年都統計為今年的第61天。 第二個錯誤來源於原文中的錯字,Michalsky在1989年的一篇文章中更正了錯誤(太陽能43(5):323)。

此代碼塊顯示違規行,註釋掉並立即跟隨糾正版本:

# leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & day >= 60
  leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & 
              day >= 60 & !(month==2 & day==60)

# oblqec <- 23.429 - 0.0000004 * time
  oblqec <- 23.439 - 0.0000004 * time

更正的sunPosition()版本

這是上面驗證的更正的代碼:

sunPosition <- function(year, month, day, hour=12, min=0, sec=0,
                    lat=46.5, long=6.5) {

    twopi <- 2 * pi
    deg2rad <- pi / 180

    # Get day of the year, e.g. Feb 1 = 32, Mar 1 = 61 on leap years
    month.days <- c(0,31,28,31,30,31,30,31,31,30,31,30)
    day <- day + cumsum(month.days)[month]
    leapdays <- year %% 4 == 0 & (year %% 400 == 0 | year %% 100 != 0) & 
                day >= 60 & !(month==2 & day==60)
    day[leapdays] <- day[leapdays] + 1

    # Get Julian date - 2400000
    hour <- hour + min / 60 + sec / 3600 # hour plus fraction
    delta <- year - 1949
    leap <- trunc(delta / 4) # former leapyears
    jd <- 32916.5 + delta * 365 + leap + day + hour / 24

    # The input to the Atronomer's almanach is the difference between
    # the Julian date and JD 2451545.0 (noon, 1 January 2000)
    time <- jd - 51545.

    # Ecliptic coordinates

    # Mean longitude
    mnlong <- 280.460 + .9856474 * time
    mnlong <- mnlong %% 360
    mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

    # Mean anomaly
    mnanom <- 357.528 + .9856003 * time
    mnanom <- mnanom %% 360
    mnanom[mnanom < 0] <- mnanom[mnanom < 0] + 360
    mnanom <- mnanom * deg2rad

    # Ecliptic longitude and obliquity of ecliptic
    eclong <- mnlong + 1.915 * sin(mnanom) + 0.020 * sin(2 * mnanom)
    eclong <- eclong %% 360
    eclong[eclong < 0] <- eclong[eclong < 0] + 360
    oblqec <- 23.439 - 0.0000004 * time
    eclong <- eclong * deg2rad
    oblqec <- oblqec * deg2rad

    # Celestial coordinates
    # Right ascension and declination
    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] + twopi
    dec <- asin(sin(oblqec) * sin(eclong))

    # Local coordinates
    # Greenwich mean sidereal time
    gmst <- 6.697375 + .0657098242 * time + hour
    gmst <- gmst %% 24
    gmst[gmst < 0] <- gmst[gmst < 0] + 24.

    # Local mean sidereal time
    lmst <- gmst + long / 15.
    lmst <- lmst %% 24.
    lmst[lmst < 0] <- lmst[lmst < 0] + 24.
    lmst <- lmst * 15. * deg2rad

    # Hour angle
    ha <- lmst - ra
    ha[ha < -pi] <- ha[ha < -pi] + twopi
    ha[ha > pi] <- ha[ha > pi] - twopi

    # Latitude to radians
    lat <- lat * deg2rad

    # Azimuth and elevation
    el <- asin(sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(ha))
    az <- asin(-cos(dec) * sin(ha) / cos(el))

    # For logic and names, see Spencer, J.W. 1989. Solar Energy. 42(4):353
    cosAzPos <- (0 <= sin(dec) - sin(el) * sin(lat))
    sinAzNeg <- (sin(az) < 0)
    az[cosAzPos & sinAzNeg] <- az[cosAzPos & sinAzNeg] + twopi
    az[!cosAzPos] <- pi - az[!cosAzPos]

    # if (0 < sin(dec) - sin(el) * sin(lat)) {
    #     if(sin(az) < 0) az <- az + twopi
    # } else {
    #     az <- pi - az
    # }


    el <- el / deg2rad
    az <- az / deg2rad
    lat <- lat / deg2rad

    return(list(elevation=el, azimuth=az))
}

參考文獻:

Michalsky,JJ 1988.用於近似太陽位置的天文年曆算法(1950-2050)。 太陽能。 40(3):227-235。

Michalsky,JJ 1989.勘誤。 太陽能。 43(5):323。

Spencer,JW 1989.關於“天文年曆的近似太陽位置算法(1950-2050)”的評論。 太陽能。 42(4):353。

Walraven,河1978年。計算太陽的位置。 太陽能。 20:393-397。


Answer #4

這是對喬希傑出答案的建議更新。

該函數的大部分開始是用於計算自2000年1月1日中午以來的天數的樣板代碼。這更好地處理了使用R的現有日期和時間函數。

我還認為,不是指定日期和時間有六個不同的變量,而是更容易(並且與其他R函數更一致)指定現有的日期對像或日期字符串+格式字符串。

這裡有兩個輔助函數

astronomers_almanac_time <- function(x)
{
  origin <- as.POSIXct("2000-01-01 12:00:00")
  as.numeric(difftime(x, origin, units = "days"))
}

hour_of_day <- function(x)
{
  x <- as.POSIXlt(x)
  with(x, hour + min / 60 + sec / 3600)
}

現在函數的開始簡化為

sunPosition <- function(when = Sys.time(), format, lat=46.5, long=6.5) {

  twopi <- 2 * pi
  deg2rad <- pi / 180

  if(is.character(when)) when <- strptime(when, format)
  time <- astronomers_almanac_time(when)
  hour <- hour_of_day(when)
  #...

另一個奇怪的是像這樣的線條

mnlong[mnlong < 0] <- mnlong[mnlong < 0] + 360

由於mnlong已經%%呼籲其值,它們都應該是非負的,所以這條線是多餘的。





azimuth