Tedward - 11 months ago 45

R Question

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

Answer Source

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`

.