CraigV
Posts: 36
Joined: Tue May 21, 2013 11:45 am

Calculation of daylight times from Earth orbital data, etc.

Wed Jun 25, 2014 7:17 am

Using Earth orbital data, location, and date, this class and sample program starts the camera at daybreak and shuts it down at nightfall. The critical times are produced by a SunriseSunset class which the remaining code uses to run the camera.
The calculation takes into account the elliptical nature of the Earth's orbit and should be good throughout the year. I run it in the background and let it scp the results to another computer. I hope it might be useful to others. I have not checked it in the Southern hemisphere so let me know if there is a problem there or any other fixes.

Correction: The sunrise_sunset class is ok, but when cleaning up the code for placement here, I introduced a bug in the sample code that used the sunrise_sunset class. That bug showed up on July 1 when the program used the output of time.asctime() to create a file name. A space was introduced where a '0' was desired. The code below is now fixed in the three lines that were defective. Two of the fixed lines start with videoFileName= and the third starts with 'stillFileName=' and all three are continued on the following line. A ':>0s' is removed and a '.replace(' ','0') is added in each case.

Code: Select all

#! /usr/bin/python3
""" This program records video from daybreak to nightfall as calculated from the earth's orbital properties.

The camera location must be provided and certain directories specified as noted in the code comments.

If put in the background, and then stopped, be sure to kill the raspivid subprocess.
"""
import math
import time
import pickle
import datetime
import subprocess

class sunrise_sunset():
    """Calculation of daybreak, sunrise, midday, sunset, and nighfall for a time and location.""" 
    def __init__(self,lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_WE):
        """Locating given as latitude degrees, latitude minutes with 'N' or 'S'."""
        """                  longidute degrees, longitude minutes with 'E' or 'W'."""
        if lat_NS=='N':
            self.latitude=lat_deg+lat_min/60.
        else:
            self.latitude=(-lat_deg-lat_min/60.)*math.pi/180
        if lon_WE=='W':
            self.longitude=lon_deg+lon_min/60.
        else:
            self.longitude=(360-lon_deg-lon_min/60.)*math.pi/180
        self.latitudeRad=self.latitude*math.pi/180
        self.longitudeRad=self.longitude*math.pi/180

    def JulianToUnixDays(t):
        """Calculate Julian day from Unix day given as a datetime.datetime object."""
        return t-2440587.5

    def julianDaysFromEpoch2000(t):
        """Calculate Julian day from datetime.datetime object."""
        if t.hour<12:
            julianOffset=-1
        else:
            julianOffset=0
        return datetime.datetime.fromordinal(t.toordinal()+julianOffset-datetime.datetime(2000,1,1,0,0,0,0,None).toordinal()).toordinal()

    def julianDaysToDatetime(Jt):
        """Calculate local datetime.datetime object from Julian day."""
        t0=sunrise_sunset.JulianToUnixDays(Jt)+datetime.date(1970,1,1).toordinal()
        t0=t0-(datetime.datetime.utcnow()-datetime.datetime.now()).seconds/86400
        t0_days=int(t0)
        t0_day_fraction=t0-t0_days
        t0_hours=int(t0_day_fraction*24)
        t0_minutes=int(t0_day_fraction*1440-t0_hours*60)
        t0_seconds=int(t0_day_fraction*86400-t0_hours*3600-t0_minutes*60)
        return datetime.datetime.combine(datetime.date.fromordinal(int(t0)),datetime.time(hour=t0_hours, minute=t0_minutes, second=t0_seconds))

    def calculate(self,now):
        """Calculate daybreak, sunrise, midday, sunset, and nightfall for current local day."""

        # Julian cycle since Jan 1, 2000
        n=sunrise_sunset.julianDaysFromEpoch2000(now)

        # Approximate Solar Noon
        Jstar=2451545.0009+self.longitude/360+n

        # Solar Mean Anomaly
        MDeg=(357.5291 + 0.98560028*(Jstar-2451545)) % 360
        MRad=MDeg*math.pi/180

        # Equation of Center
        CDeg = 1.9148*math.sin(MRad) + 0.02*math.sin(2*MRad) + 0.0003*math.sin(3*MRad)

        # Ecliptic Longitude
        LambdaDeg=(MDeg + 102.9372 + CDeg + 180) % 360
        LambdaRad=LambdaDeg*math.pi/180

        # Solar Transit
        Jtransit=Jstar + 0.0053*math.sin(MRad) - 0.0069*math.sin(2*LambdaRad)

        # Earth Axis Tilt
        tiltDeg=23.45
        tiltRad=tiltDeg*math.pi/180

        # Declination of the Sun
        sinDelta=math.sin(LambdaRad)*math.sin(tiltRad)
        cosDelta=math.cos(math.asin(sinDelta))
        deltaRad=math.asin(sinDelta)
        deltaDeg=deltaRad*180/math.pi

        # Hour Angle
        omega0Rad=math.acos((math.sin(-0.83*math.pi/180) - math.sin(self.latitudeRad)*sinDelta)/(math.cos(self.latitudeRad)*cosDelta))
        omega0Deg=omega0Rad*180/math.pi

        # Sunet and Sunrise
        Jset=2451545.0009 + (omega0Deg+self.longitude)/360 + n + 0.0053*math.sin(MRad)-0.0069*math.sin(2*LambdaRad)
        Jrise=Jtransit-(Jset-Jtransit)

        self.offset=abs(int((Jtransit-Jset)*75))  # This is the difference in minutes between daybreak and sunrise
        afterJset=Jset+self.offset/1440.          #   or between sunset and nightfall.  The value of 75 was
        beforeJrise=Jrise-self.offset/1440.       #   determined by when there was too little light for the camera.  

        self.daybreak=sunrise_sunset.julianDaysToDatetime(beforeJrise)
        self.sunrise=sunrise_sunset.julianDaysToDatetime(Jrise)
        self.transit=sunrise_sunset.julianDaysToDatetime(Jtransit)
        self.sunset=sunrise_sunset.julianDaysToDatetime(Jset)
        self.nightfall=sunrise_sunset.julianDaysToDatetime(afterJset)

# These directories need to be available and need to be manually cleared as appropriate.
basicVideoOptions='--nopreview --metering matrix --inline --awb sun --bitrate 0 --exposure auto'
basicImageOptions='--exposure auto --awb sun'
localVideoDir='/home/pi/Camera/Videos/'
localImageDir='/home/pi/Camera/Images/'
# A remote computer where ssh is accepted by key exchange for remote storage.
remoteVideoDir='marren@marren:/media/Maxtor-1/Camera/Videos/'
remoteImageDir='marren@marren:/media/Maxtor-1/Camera/Images/'

timingCorrectionFactor=0.9986  # Correction factor for video time.  Pi video clock appears to run 0.14% too fast.

maxSessionLength=7200 # seconds

# This log file needs to be manually cleared as appropriate.
logFileName='CameraLog'
def logNote(s):
    logFile=open(logFileName,'a')
    logFile.write(s)
    logFile.close()

# The following needs to be set to the correct location.
# House front gate 37 degrees 14.321' N
#                 119 degrees 55.037' W
s=sunrise_sunset(37, 14.321, 'N', 119, 55.037, 'W')
logNote('Calculations are for {:8.5f} north latitude\n'
        '                and {:9.5f} west longitude.\n'.format(s.latitude, s.longitude))

# Repeat taking pictures forever but note that files must be deleted manually
#    every few days to prevent SD disk becoming filled.
# If the program is run without a terminal (in the background) and killed when the video is running, a reboot will be necessary.
# If it is killed after nightfall and before daybreak, no reboot should be needed.
while True:
    remoteStorageFailed=False
    currentDateTime = datetime.datetime.now()
    s.calculate(currentDateTime)
    logNote('currentDateTime={}\n'.format(currentDateTime)+\
            '   s.offset={} minutes\n'.format(s.offset))

    if currentDateTime > s.nightfall:
        tomorrowDateTime=datetime.datetime.combine(datetime.datetime.fromordinal(datetime.date.today().toordinal()+1),datetime.datetime.now().time())
        s.calculate(tomorrowDateTime)

    logNote(' s.daybreak={}\n'.format(s.daybreak)+\
            '  s.sunrise={}\n'.format(s.sunrise)+\
            '  s.transit={}\n'.format(s.transit)+\
            '   s.sunset={}\n'.format(s.sunset)+\
            's.nightfall={}\n'.format(s.nightfall))

    if s.daybreak < currentDateTime < s.transit:
        am_sessions_length=s.transit - currentDateTime
        pm_sessions_length=s.nightfall-s.transit
    elif s.transit < currentDateTime < s.nightfall:
        am_sessions_length=datetime.timedelta()
        pm_sessions_length=s.nightfall - currentDateTime
    else:
        am_sessions_length=s.transit-s.daybreak
        pm_sessions_length=s.nightfall-s.transit

    if am_sessions_length.seconds > 0:
        am_sessions=int(am_sessions_length.seconds/maxSessionLength) + 1
        am_session_length=am_sessions_length/am_sessions
        logNote('\n  Morning: {} sessions of length {} ms for a total of {}\n'\
                .format(am_sessions,am_session_length.seconds*1000,am_sessions_length))

    if pm_sessions_length.seconds > 0:
        pm_sessions=int(pm_sessions_length.seconds/maxSessionLength) + 1
        pm_session_length=pm_sessions_length/pm_sessions
        logNote('Afternoon: {} sessions of length {} ms for a total of {}\n'\
                .format(pm_sessions,pm_session_length.seconds*1000,pm_sessions_length))

    # Morning session and still image at midday
    if currentDateTime < s.transit:
        if currentDateTime < s.daybreak:
            waitDuration=(s.daybreak-currentDateTime).seconds
            logNote('Waiting {} seconds...\n'.format(waitDuration))
            time.sleep(waitDuration)

        for session in range(am_sessions):
            if session==0 or session==1:
                qp=25
                width=1920
                height=1080
                fps=2
            else:
                qp=30
                width=1296
                height=730
                fps=2
            t=time.asctime()
            videoFileName='{}-{}-{}-{}-gp{}-{}-{}x{}x{}'\
                 .format(t[4:7],t[8:10].replace(' ','0'),t[20:],t[11:13]+t[14:16]+t[17:19],qp,str(am_session_length.seconds),width,height,fps)
            videoCommand='/usr/bin/raspivid {} --qp {} --width {} --height {} --timeout {} --framerate {} --output {}.h264'\
                 .format(basicVideoOptions,qp,width,height,int(am_session_length.seconds*timingCorrectionFactor*1000),fps,localVideoDir+videoFileName)
            logNote(videoCommand+'\n')
            # Run the next video session and wait for it to finish.
            p = subprocess.call(videoCommand.split())

            if not remoteStorageFailed:
                try:
                    transferCommand='/usr/bin/scp -q {}.h264 {}.h264'\
                                .format(localVideoDir+videoFileName, remoteVideoDir+videoFileName)
                    # Popen is used here to start the transfer and return without waiting.
                    subprocess.Popen(transferCommand.split())
                    logNote('{} transfer started.\n'.format(videoFileName))
                except:
                    remoteStorageFailed=True 

        # Take a still at midday
        t=time.asctime()
        stillFileName='{}-{}-{}-{}-2592x1944'\
             .format(t[4:7],t[8:10].replace(' ','0'),t[20:],t[11:13]+t[14:16]+t[17:19])
        imageCommand='/usr/bin/raspistill {} --width 2592 --height 1944 --timeout 5000 --output {}.jpg'\
             .format(basicImageOptions, localImageDir+stillFileName)
        logNote(imageCommand+'\n')
        p = subprocess.call(imageCommand.split())

        if not remoteStorageFailed:
            try:
                transferCommand="/usr/bin/scp -q {}.jpg {}.jpg"\
                    .format(localImageDir+stillFileName, remoteImageDir+stillFileName)
                # Popen is used here to start the transfer and return without waiting.
                subprocess.Popen(transferCommand.split())
                logNote('{} transfer started.\n'.format(videoFileName))
            except:
                remoteStorageFailed=True 
                logNote('Remote storage attempt failed..\n')

    # Afternoon session
    remoteStorageFailed=False
    for session in range(pm_sessions):
        if session==pm_sessions-1:
            qp=25
            width=1920
            height=1080
            fps=2
        else:
            qp=30
            width=1296
            height=730
            fps=2
        t=time.asctime()
        videoFileName='{}-{}-{}-{}-gp{}-{}-{}x{}x{}'\
             .format(t[4:7],t[8:10].replace(' ','0'),t[20:],t[11:13]+t[14:16]+t[17:19],qp,str(pm_session_length.seconds),width,height,fps)
        videoCommand='/usr/bin/raspivid {} --qp {} --width {} --height {} --timeout {} --framerate {} --output {}.h264'\
             .format(basicVideoOptions,qp,width,height,int(pm_session_length.seconds*timingCorrectionFactor*1000),fps,localVideoDir+videoFileName)
        logNote(videoCommand+'\n')
        # Run the next video session and wait for it to finish.
        p = subprocess.call(videoCommand.split())

        if not remoteStorageFailed:
            try:
                transferCommand="/usr/bin/scp -q {}.h264 {}.h264"\
                    .format(localVideoDir+videoFileName, remoteVideoDir+videoFileName)
                # Popen is used here to start the transfer and return without waiting.
                subprocess.Popen(transferCommand.split())
                logNote('{} transfer started.\n'.format(videoFileName))
            except:
                remoteStorageFailed=True 
                logNote('Remote storage attempt failed..\n')

    # Show files
    output=subprocess.check_output(['/bin/ls','-l','Videos/'])
    logNote(output.decode())

    output=subprocess.check_output(['/bin/ls','-l','Images/'])
    logNote(output.decode())

    # Wait 10 minutes before preparing for the next day.
    time.sleep(600)
Last edited by CraigV on Tue Jul 01, 2014 9:04 pm, edited 1 time in total.

User avatar
DougieLawson
Posts: 42142
Joined: Sun Jun 16, 2013 11:19 pm
Location: A small cave in deepest darkest Basingstoke, UK

Re: Calculation of daylight times from Earth orbital data, e

Wed Jun 25, 2014 7:23 am

http://www.raspberrypi.org/forums/viewt ... 91&t=71010

Read to the bottom of that thread, there's a program called sunwait that can be set-up in crontab to run at sunrise and sunset (or in my case at civil twilight).

https://www.risacher.org/sunwait/




How's that for my 5000th forum post?
Languages using left-hand whitespace for syntax are ridiculous

DMs sent on https://twitter.com/DougieLawson or LinkedIn will be answered next month.
Fake doctors - are all on my foes list.

The use of crystal balls and mind reading is prohibited.

tenochtitlanuk
Posts: 156
Joined: Fri Jul 06, 2012 8:51 pm
Location: Taunton, Somerset, UK

Re: Calculation of daylight times from Earth orbital data, e

Wed Jun 25, 2014 8:06 am

Just a word of caution. This code implements the algorithm shown on Wikipedia. Be aware that for some sites it fails.
There are subtle errors related to altitude and atmospheric conditions - rarely worth worrying about - but much greater errors affecting high latitudes. If you are inside the arctic circle then the 'Wikipedia' algorithm is useless. ( On the first day of Arctic 'Summer' there's a sunrise but no sunset, for example.)
We DO have Pi owners in strange places!!

CraigV
Posts: 36
Joined: Tue May 21, 2013 11:45 am

Re: Calculation of daylight times from Earth orbital data, e

Wed Jun 25, 2014 8:31 am

Yes. I forgot to mention the problem near the poles. I lump the altitude problem together with considerations of local terrain. Whatever its weaknesses, it has worked well for me since October, 2013, through both solstices.

Viila
Posts: 3
Joined: Mon Oct 21, 2013 12:28 pm

Re: Calculation of daylight times from Earth orbital data, e

Wed Jun 25, 2014 7:09 pm

You might want to consider also something called PyEphem (http://rhodesmill.org/pyephem/), which should be able to provide you with all the basic astronomical calculations that a normal healthy citizen needs ;)

Return to “Camera board”