In [1]:

    

%run basics
%matplotlib
import ephem
import pysolar
import pytz


    

    




Using matplotlib backend: Qt4Agg

In [2]:

    

def solar_zenith_im(ldt,latitude,longitude,time_zone):
    # get a pytz object
    tz = pytz.timezone(time_zone)
    # use the same "current" date as Ian
    cdt = datetime.datetime(2013,6,1,0,0,0)
    # get the offset to UTC time as a timedelta object
    utcoffset = tz.utcoffset(cdt)
    # get the offset as hours
    GMT_zone = utcoffset.seconds/float(3600)
    # get the day of the year from local time
    DOY = numpy.array([dt.timetuple().tm_yday for dt in ldt])
    # For each day calculate equation of time correction, solar noon and declination
    array_EqofTime=0.17*numpy.sin(4*numpy.pi*(DOY-80)/373)-0.129*numpy.sin(2*numpy.pi*(DOY-8)/355) # DiLaura (1984)
    array_solar_noon=12+(GMT_zone*15.0-longitude)/360*24-array_EqofTime # Me
    array_decl=numpy.radians(23.4)*numpy.sin((DOY+284)/365.0*2*numpy.pi) # Oke (1987)
    # Create an hour angle array for each minute of day and each day of year
    array_h=numpy.tile(numpy.linspace(0,1439.0/1440*24,num=1440),(len(DOY),1))
    array_h=abs(numpy.radians((array_solar_noon.reshape(len(DOY),1)-array_h)*15))
    # Duplicate declination array for each time of day
    array_decl=numpy.tile(array_decl,(1440,1)).T
    # Calculate zenith angles
    array_z=numpy.arccos(numpy.sin(numpy.radians(latitude))*numpy.sin(array_decl)+numpy.cos(numpy.radians(latitude))*numpy.cos(array_decl)*numpy.cos(array_h))
    return array_z,array_h


    

In [3]:

    

def solar_zenith_pri(ldt,latitude,longitude,time_zone):
    # get a pytz object
    tz = pytz.timezone(time_zone)
    # use the same "current" date as Ian
    cdt = datetime.datetime(2013,6,1,0,0,0)
    # get the offset to UTC time as a timedelta object
    utcoffset = tz.utcoffset(cdt)
    # get the offset as hours
    GMT_zone = utcoffset.seconds/float(3600)
    # get the day of the year from local time
    DOY = numpy.array([dt.timetuple().tm_yday for dt in ldt])
    # For each day calculate equation of time correction, solar noon and declination
    eqn_of_time = 0.17*numpy.sin(4*numpy.pi*(DOY-80)/373)-0.129*numpy.sin(2*numpy.pi*(DOY-8)/355)
    solar_noon = 12+(GMT_zone*15.0-site_longitude)/360*24-eqn_of_time
    decln = numpy.radians(23.4)*numpy.sin((DOY+284)/365.0*2*numpy.pi)
    Hdh = numpy.array([dt.hour+dt.minute/float(60) for dt in ldt])
    Hdh = abs(numpy.radians((solar_noon-Hdh)*15))
    solar_zenith = numpy.arccos(numpy.sin(numpy.radians(site_latitude))
                                *numpy.sin(decln)
                                +numpy.cos(numpy.radians(site_latitude))
                                *numpy.cos(decln)*numpy.cos(Hdh))
    return solar_zenith,Hdh


    

In [4]:

    

def solar_altitude_from_zenith(solar_zenith):
    return numpy.pi/2-solar_zenith


    

In [5]:

    

nc_name = "../Sites/HowardSprings/Data/Portal/HowardSprings_2014_L3.nc"


    

In [6]:

    

ds = qcio.nc_read_series(nc_name)
site_name = ds.globalattributes["site_name"]
site_timezone = ds.globalattributes["time_zone"]
site_latitude = float(ds.globalattributes["latitude"])
site_longitude = float(ds.globalattributes["longitude"])
ldt = ds.series["DateTime"]["Data"]


    

In [7]:

    

tz = pytz.timezone(site_timezone)
cdt = datetime.datetime(2013,6,1,0,0,0)
utcoffset = tz.utcoffset(cdt)
GMT_zone = utcoffset.seconds/float(3600)
DOY = numpy.array([dt.timetuple().tm_yday for dt in ldt])


    

In [8]:

    

sz_pri,Hdh = solar_zenith_pri(ldt,site_latitude,site_longitude,site_timezone)
sz_im,array_h = solar_zenith_im(ldt,site_latitude,site_longitude,site_timezone)


    

In [ ]:

    

fig=plt.figure()
plt.plot(Hdh[0:47],sz_pri[0:47])
plt.show()


    

In [ ]:

    

fig=plt.figure()
plt.plot(array_h[0,:],sz_im[0,:])
plt.show()


    

In [9]:

    

sa_pri = solar_altitude_from_zenith(sz_pri)
sa_im = solar_altitude_from_zenith(sz_im)


    

In [10]:

    

# get the UTC time from the local time
ldt_UTC = qcutils.get_UTCfromlocaltime(ds)
# get the solar altitude
sa_pysolar = [pysolar.GetAltitude(site_latitude,site_longitude,dt) for dt in ldt_UTC]


    

In [26]:

    

obs=ephem.Observer()
obs.lat = str(site_latitude)
obs.lon = str(site_longitude)
sa_ephem = numpy.zeros(len(ldt_UTC))
for i in range(0,len(ldt_UTC)):
    obs.date = ldt_UTC[i]
    sun = ephem.Sun(obs)
    sun.compute(obs)
    sa_ephem[i] = numpy.rad2deg(float(sun.alt))


    

In [29]:

    

for i in range(0,48):
    print numpy.rad2deg(sa_pri[i]),sa_pysolar[i],numpy.rad2deg(sa_pri[i])-sa_pysolar[i],sa_ephem[i],sa_ephem[i]-sa_pysolar[i]


    

    




-54.2334862534 -54.1771478237 -0.0563384297386 -54.1670850053 0.0100628183951
-54.442946223 -54.3569084923 -0.0860377306191 -54.3463535653 0.0105549270378
-53.159151003 -53.0504466984 -0.108704304623 -53.0387172638 0.011729434557
-50.5278488378 -50.4049790029 -0.122869834877 -50.3914657151 0.0135132877972
-46.7951760718 -46.6658239231 -0.129352148623 -46.6499461444 0.0158777787714
-42.2194642768 -42.0895306276 -0.129933649139 -42.0706219097 0.0189087179598
-37.0198424912 -36.8936452987 -0.126197192568 -36.8708296304 0.022815668267
-31.3634298969 -31.2444096147 -0.119020282117 -31.2163550732 0.0280545415681
-25.3716658899 -25.2633117335 -0.108354156435 -25.2277664448 0.0355452887133
-19.1314293597 -19.0386825765 -0.0927467832255 -18.9914161193 0.0472664571867
-12.7050472593 -12.6367514101 -0.0682958491953 -12.5693119692 0.0674394409335
-6.13784016062 -6.06042895103 -0.0774112095906 -5.83412603476 0.226302916269
0.536544586259 1.03974729543 -0.503202709174 1.05181690223 0.0120696067942
7.2924221445 7.5237331743 -0.231311029797 7.52707080971 0.00333763541095
14.1094445637 14.2855904026 -0.176145838891 14.2877858568 0.00219545424543
20.9706408411 21.1216652807 -0.151024439605 21.124137565 0.00247228428687
27.860784047 27.9958297667 -0.135045719701 27.9981492287 0.00231946203337
34.7647513642 34.8876037147 -0.122852350455 34.8899040164 0.00230030168991
41.6654169821 41.7777912022 -0.11237422008 41.780086153 0.00229495074007
48.5402076763 48.6426994338 -0.102491757484 48.644972684 0.00227325017563
55.3542753205 55.4465650076 -0.092289687141 55.4487768332 0.00221182560328
62.0446589503 62.1252236895 -0.0805647391288 62.1273050678 0.00208137832811
68.4775773158 68.5425796351 -0.0650023192803 68.5444043968 0.00182476172661
74.3148548413 74.3549346368 -0.0400797954864 74.3562055314 0.00127089458138
78.5854044944 78.5798188026 0.00558569184922 78.5798510801 3.22775789243e-05
79.2198470522 79.1568505453 0.0629965068341 79.1550759517 -0.00177459367085
75.7240511664 75.6384973164 0.0855538499895 75.6357936611 -0.00270365538422
70.1777064027 70.0924164512 0.0852899515207 70.0894751493 -0.00294130185004
63.8607808865 63.7811302851 0.0796506014167 63.7781412996 -0.00298898548814
57.2233659299 57.1505050502 0.0728608797414 57.1475029364 -0.00300211378272
50.4353889765 50.3697450368 0.0656439397177 50.3667404303 -0.00300460650509
43.5731389626 43.5152160083 0.0579229543632 43.5122016359 -0.00301437235711
36.6770250304 36.6277197284 0.0493053019991 36.6247003486 -0.0030193798539
29.772105705 29.7330780599 0.0390276451876 29.7300717588 -0.00300630102812
22.8765798458 22.8510794709 0.0255003749329 22.848184887 -0.00289458390014
16.0058785272 16.0011786369 0.00469989035986 15.9988814326 -0.00229720428761
9.17502609502 9.21247265693 -0.0374465619119 9.20992526081 -0.00254739611367
2.40035329095 2.59660906971 -0.196255778759 2.60139942159 0.00479035187591
-4.29890363798 -4.24675350942 -0.0521501285541 -2.8832806883 1.36347282112
-10.8987656092 -11.0149285687 0.116162959462 -10.9458810176 0.0690475510698
-17.3681069078 -17.4575477958 0.0894408880172 -17.4111903831 0.0463574126805
-23.6652936502 -23.7373261753 0.0720325251128 -23.7041544391 0.0331717362331
-29.7333599061 -29.7928201078 0.0594602016991 -29.7677402521 0.0250798557118
-35.4931273976 -35.5423541815 0.049226783944 -35.5226903071 0.0196638743831
-40.8338490523 -40.87401299 0.0401639376496 -40.8581608792 0.0158521108041
-45.6020749126 -45.6337623739 0.0316874613131 -45.6206161459 0.013146228021
-49.5929745554 -49.6165332292 0.0235586738885 -49.6052187341 0.0113144951641
-52.5972258946 -52.5717427112 -0.0254831833378 -52.5614612106 0.0102815006053

In [22]:

    

print ldt[0],ldt_UTC[0],ldt_UTC[0]+datetime.timedelta(hours=9.5)


    

    




 2014-01-01 00:30:00 2013-12-31 15:00:00+00:00 2014-01-01 00:30:00+00:00

In [ ]: