4

I have UTC time (hours, minutes, seconds), longitude(deg E), latitude (deg N) and date. Can anyone provide me with a code to calculate solar zenith angle in Python 2.7?

alecxe
  • 462,703
  • 120
  • 1,088
  • 1,195
hima Fernando
  • 137
  • 2
  • 8

2 Answers2

11

It is an interesting problem and I think I have a good answer - well, at least starting points.

Check out the awesome astropy package. I believe you need to use the coordinates module.

Something along these lines:

import astropy.coordinates as coord
from astropy.time import Time
import astropy.units as u


loc = coord.EarthLocation(lon=0.1 * u.deg,
                          lat=51.5 * u.deg)
now = Time.now()

altaz = coord.AltAz(location=loc, obstime=now)
sun = coord.get_sun(now)

print(sun.transform_to(altaz).alt)

Here, we are getting the angle of the sun above the horizon for the 0.1 degrees longitude and 51.5 latitude location at the current time.

FYI, .zen would give you the zenith angle.

Community
  • 1
  • 1
alecxe
  • 462,703
  • 120
  • 1,088
  • 1,195
  • Why do they use their own `Time` class instead of the standard `datetime`? – Mark Ransom Dec 05 '17 at 22:31
  • @MarkRansom I think the `AltAz` expects a `Time()` instance, but not sure - still exploring the powers of the package. Thanks. – alecxe Dec 05 '17 at 22:32
  • Thank you for your prompt response Alecxe. I hope this will be helpful for me. Thanks again – hima Fernando Dec 05 '17 at 23:52
  • @MarkRansom - `Time` is used instead of `datetime` because `datetime` is missing many astronomy-critical features. E.g., dynamic range of a nanosecond over Gigayear spans, conversion to/from standard astronomy forms like JD or MJD, an understanding of the UT1 to UTC conversion, numpy support, etc. – eteq Dec 06 '17 at 15:09
  • Also if for whatever reason you want to use `datetime`, you can pass it in as the obstime argument and it will auto-convert to a `Time` internally. – eteq Dec 06 '17 at 15:24
  • @eteq after I asked I started speculating it was something like that, thanks. – Mark Ransom Dec 06 '17 at 15:40
7

@alecxe's answer is great, but I thought I'd add a slight modification that's a bit closer to what the original question is asking (the zenithal angle at a specific time)

from astropy.coordinates import get_sun, AltAz, EarthLocation
from astropy.time import Time

sun_time = Time('2017-12-6 17:00') #UTC time
loc = EarthLocation.of_address('Baltimore, MD')  # anything the google geocoding API resolves
altaz = AltAz(obstime=sun_time, location=loc)

zen_ang = get_sun(sun_time).transform_to(altaz).zen

zen_ang is then an Angle object - see more about those in the docs, but basically they end up working like numpy scalars with associated units of "degrees".

eteq
  • 960
  • 7
  • 6