24

I want to classify timestamp data types in a PostgreSQL table with regards to whether they can be considered "at day" or "at night". In other words I want to be able to calculate sunrise and sunset times accurately, given a particular GPS position.

I know plpgsql and plpython.

1
  • 1
    Going by accepted answer it seems that title can be changed from "within PostgreSQL" to "in Python". Is this a correct interpretation? Dec 1, 2018 at 8:29

6 Answers 6

29

Take a look at these links:

19

I know this is yonks old, but I thought I'd share since I found no quick solution. This uses the Sun class (see below), which I constructed by following this link.

from Sun import Sun

coords = {'longitude' : 145, 'latitude' : -38 }

sun = Sun()

# Sunrise time UTC (decimal, 24 hour format)
print sun.getSunriseTime( coords )['decimal']

# Sunset time UTC (decimal, 24 hour format)
print sun.getSunsetTime( coords )['decimal']

It seems to be accurate to within a few minutes, at least where I live. For greater accuracy, the zenith param in the calcSunTime() method could use fine tuning. See the above link for more info.

# save this as Sun.py

import math
import datetime

class Sun:

    def getSunriseTime( self, coords ):
        return self.calcSunTime( coords, True )

    def getSunsetTime( self, coords ):
        return self.calcSunTime( coords, False )

    def getCurrentUTC( self ):
        now = datetime.datetime.now()
        return [ now.day, now.month, now.year ]
    
    def calcSunTime( self, coords, isRiseTime, zenith = 90.8 ):
    
        # isRiseTime == False, returns sunsetTime
    
        day, month, year = self.getCurrentUTC()
    
        longitude = coords['longitude']
        latitude = coords['latitude']

        TO_RAD = math.pi/180
    
        #1. first calculate the day of the year
        N1 = math.floor(275 * month / 9)
        N2 = math.floor((month + 9) / 12)
        N3 = (1 + math.floor((year - 4 * math.floor(year / 4) + 2) / 3))
        N = N1 - (N2 * N3) + day - 30

        #2. convert the longitude to hour value and calculate an approximate time
        lngHour = longitude / 15
        
        if isRiseTime:
            t = N + ((6 - lngHour) / 24)
        else: #sunset
            t = N + ((18 - lngHour) / 24)
    
        #3. calculate the Sun's mean anomaly
        M = (0.9856 * t) - 3.289
    
        #4. calculate the Sun's true longitude
        L = M + (1.916 * math.sin(TO_RAD*M)) + (0.020 * math.sin(TO_RAD * 2 * M)) + 282.634
        L = self.forceRange( L, 360 ) #NOTE: L adjusted into the range [0,360)
    
        #5a. calculate the Sun's right ascension
    
        RA = (1/TO_RAD) * math.atan(0.91764 * math.tan(TO_RAD*L))
        RA = self.forceRange( RA, 360 ) #NOTE: RA adjusted into the range [0,360)
    
        #5b. right ascension value needs to be in the same quadrant as L
        Lquadrant  = (math.floor( L/90)) * 90
        RAquadrant = (math.floor(RA/90)) * 90
        RA = RA + (Lquadrant - RAquadrant)
    
        #5c. right ascension value needs to be converted into hours
        RA = RA / 15
    
        #6. calculate the Sun's declination
        sinDec = 0.39782 * math.sin(TO_RAD*L)
        cosDec = math.cos(math.asin(sinDec))
    
        #7a. calculate the Sun's local hour angle
        cosH = (math.cos(TO_RAD*zenith) - (sinDec * math.sin(TO_RAD*latitude))) / (cosDec * math.cos(TO_RAD*latitude))
    
        if cosH > 1:
            return {'status': False, 'msg': 'the sun never rises on this location (on the specified date)'}
    
        if cosH < -1:
            return {'status': False, 'msg': 'the sun never sets on this location (on the specified date)'}
    
        #7b. finish calculating H and convert into hours
        
        if isRiseTime:
            H = 360 - (1/TO_RAD) * math.acos(cosH)
        else: #setting
            H = (1/TO_RAD) * math.acos(cosH)
        
        H = H / 15
    
        #8. calculate local mean time of rising/setting
        T = H + RA - (0.06571 * t) - 6.622
    
        #9. adjust back to UTC
        UT = T - lngHour
        UT = self.forceRange( UT, 24) # UTC time in decimal format (e.g. 23.23)
    
        #10. Return
        hr = self.forceRange(int(UT), 24)
        min = round((UT - int(UT))*60,0)
    
        return {
            'status': True,
            'decimal': UT,
            'hr': hr,
            'min': min 
        }

    def forceRange( self, v, max ):
        # force v to be >= 0 and < max
        if v < 0:
            return v + max
        elif v >= max:
            return v - max
    
        return v
3
12

Use Astral (current version 1.6). The first example in the documentation shows the calculation of sunrise and sunset for a given location. A simpler example with custom latitude and longitude would be:

from datetime import date
import astral
loc = astral.Location(('Bern', 'Switzerland', 46.95, 7.47, 'Europe/Zurich', 510))
for event, time in loc.sun(date.today()).items():
    print(event, 'at', time)

Gives:

noon at 2018-03-12 12:39:59+01:00
sunset at 2018-03-12 18:30:11+01:00
sunrise at 2018-03-12 06:49:47+01:00
dusk at 2018-03-12 20:11:39+01:00
dawn at 2018-03-12 05:08:18+01:00

Then you can maybe use this as a starting point for writing your own postgres (or postgis) functions using plpython instead of plr.

2
  • for any location you can call Astral().sun_utc(date.today(), latitude, longitude) or variations for local time, different twilight periods, etc
    – rcoup
    Dec 21, 2018 at 9:34
  • 1
    astral API changed quite a bit since 1.6 but still a nice tool - check out the docs linked above for example usage. Dec 31, 2021 at 19:17
6

So I know this is an ages old question, but I keep needing to calculate this actually within Postgres. I've therefore ported oortCloud's answer over to PL/pgSQL. Perhaps it will be useful to someone:

CREATE OR REPLACE FUNCTION FORCE_RANGE(
    v DOUBLE PRECISION,
    max DOUBLE PRECISION
) RETURNS DOUBLE PRECISION AS $$
BEGIN
    IF v < 0 THEN
        RETURN v + max;
    ELSEIF v >= max THEN
        return v - max;
    END IF;

    return v;
END; $$
LANGUAGE plpgsql IMMUTABLE;


CREATE OR REPLACE FUNCTION RISE_SET_TIME(
    latitude DOUBLE PRECISION,
    longitude DOUBLE PRECISION,
    isRiseTime BOOL,
    as_of TIMESTAMPTZ,
    zenith DOUBLE PRECISION DEFAULT 90.8
)
RETURNS TIMESTAMPTZ AS $$
    DECLARE as_of_utc TIMESTAMPTZ;
    DECLARE as_of_year INT;
    DECLARE as_of_month INT;
    DECLARE as_of_day INT;

    DECLARE N1 INT;
    DECLARE N2 INT;
    DECLARE N3 INT;
    DECLARE N INT;

    DECLARE longitude_hour DOUBLE PRECISION;
    DECLARE M DOUBLE PRECISION;

    DECLARE t DOUBLE PRECISION;
    DECLARE L DOUBLE PRECISION;
    DECLARE RA DOUBLE PRECISION;

    DECLARE Lquadrant INT;
    DECLARE RAquadrant INT;
    DECLARE sinDec DOUBLE PRECISION;
    DECLARE cosDec DOUBLE PRECISION;
    DECLARE cosH DOUBLE PRECISION;
    DECLARE H DOUBLE PRECISION;
    DECLARE UT DOUBLE PRECISION;

    DECLARE hr INT;
    DECLARE min INT;
BEGIN
    as_of_utc = as_of at time zone 'utc';
    as_of_year = EXTRACT(YEAR FROM as_of_utc);
    as_of_month = EXTRACT(MONTH FROM as_of_utc);
    as_of_day = EXTRACT(DAY FROM as_of_utc);

    -- 1. first calculate the day of the year
    N1 = FLOOR(275.0 * as_of_month / 9.0);
    N2 = FLOOR((as_of_month + 9) / 12.0);
    N3 = (1 + FLOOR((as_of_year - 4 * FLOOR(as_of_year / 4.0) + 2) / 3.0));
    N = N1 - (N2 * N3) + as_of_day - 30;

    -- 2. convert the longitude to hour value and calculate an approximate time
    longitude_hour = longitude / 15.0;

    IF isRiseTime THEN
        t = N + ((6 - longitude_hour) / 24.);
    ELSE
        t = N + ((18 - longitude_hour) / 24.);
    END IF;

    -- 3. calculate the Sun's mean anomaly
    M = (0.9856 * t) - 3.289;

    -- 4. calculate the Sun's true longitude
    L = M + (1.916 * SIN(RADIANS(M))) + (0.020 * SIN(RADIANS(2 * M))) + 282.634;
    -- NOTE: L adjusted into the range [0,360)
    L = FORCE_RANGE(L, 360.0);

    -- 5a. calculate the Sun's right ascension
    RA = (1/RADIANS(1)) * ATAN(0.91764 * TAN(RADIANS(L)));
    RA = FORCE_RANGE( RA, 360 );  -- NOTE: RA adjusted into the range [0,360);

    -- 5b. right ascension value needs to be in the same quadrant as L
    Lquadrant = FLOOR(L/90.) * 90;
    RAquadrant = FLOOR(RA/90.) * 90;
    RA = RA + (Lquadrant - RAquadrant);

    -- 5c. right ascension value needs to be converted into hours
    RA = RA / 15.0;

    -- 6. calculate the Sun's declination
    sinDec = 0.39782 * SIN(RADIANS(L));
    cosDec = COS(ASIN(sinDec));

    -- 7a. calculate the Sun's local hour angle
    cosH = (COS(RADIANS(zenith)) - (sinDec * SIN(RADIANS(latitude)))) / (cosDec * COS(RADIANS(latitude)));

    IF cosH > 1 THEN
        RAISE NOTICE 'The sun never rises on this location on the specified date';
        RETURN NULL;
    END IF;

    IF cosH < -1 THEN
        RAISE NOTICE 'The sun never sets on this location on the specified date';
        RETURN NULL;
    END IF;

    -- 7b. finish calculating H and convert into hours
    IF isRiseTime THEN
        H = 360 - (1/RADIANS(1)) * ACOS(cosH);
    ELSE
        H = (1/RADIANS(1)) * ACOS(cosH);
    END IF;

    H = H / 15.0;

    -- calculate local mean time of rising/setting
    T = H + RA - (0.06571 * t) - 6.622;

    -- 9. adjust back to UTC
    UT = T - longitude_hour;
    UT = FORCE_RANGE( UT, 24);  -- UTC time in decimal format (e.g. 23.23)

    -- 10. Return
    hr = FORCE_RANGE(UT::INT, 24);
    min = ROUND((UT - UT::INT) * 60);

--     Enable for debugging purposes:
--     RAISE NOTICE 'as_of_utc: %', as_of_utc;
--     RAISE NOTICE 'as_of_year: %', as_of_year;
--     RAISE NOTICE 'as_of_month: %', as_of_month;
--     RAISE NOTICE 'as_of_day: %', as_of_day;
--     RAISE NOTICE 'N1: %', N1;
--     RAISE NOTICE 'N2: %', N2;
--     RAISE NOTICE 'N3: %', N3;
--     RAISE NOTICE 'N: %', N;
--     RAISE NOTICE 'longitude_hour: %', longitude_hour;
--     RAISE NOTICE 'M: %', M;
--     RAISE NOTICE 't: %', t;
--     RAISE NOTICE 'L: %', L;
--     RAISE NOTICE 'RA: %', RA;
--     RAISE NOTICE 'Lquadrant: %', Lquadrant;
--     RAISE NOTICE 'RAquadrant: %', RAquadrant;
--     RAISE NOTICE 'sinDec: %', sinDec;
--     RAISE NOTICE 'cosDec: %', cosDec;
--     RAISE NOTICE 'cosH: %', cosH;
--     RAISE NOTICE 'H: %', H;
--     RAISE NOTICE 'UT: %', UT;
--     RAISE NOTICE 'hr: %', hr;
--     RAISE NOTICE 'min: %', min;

    return as_of_utc::DATE + (INTERVAL '1 hour' * hr) + (INTERVAL '1 minute' * min);
END; $$
LANGUAGE plpgsql IMMUTABLE;

Example use:

SELECT
       RISE_SET_TIME(39.399872, -8.224454, TRUE, NOW()) AS rise,
       RISE_SET_TIME(39.399872, -8.224454, FALSE, NOW()) AS set
;
2

I use this to calculate the Sunrise, Sunset, Dawn and Dusk times.

Just need to replace the X's with your Coordinates. Note that the times returned are UTC so you need to add your particular time zone hours if needed.

request = Request('http://api.sunrise-sunset.org/json?lat=-XX.XXXXX&lng=XX.XXXXX&formatted=0')
response = urlopen(request)
timestring = response.read()

utcsunrise = timestring[34:39]
utcsunset = timestring[71:76]
utcmorning = timestring[182:187]
utcnight = timestring[231:236]
1
  • 5
    Please at least include your imports. Feb 6, 2018 at 15:11
-3

I'm using Sun.py. Today I got the value of minutes = 60, UT = 12.9979740551

class Sun:

#10. Return

        min = round((UT - int(UT))*60,0)

#add this after min calculate

        if min == 60:
            hr += 1
            min = 0
1
  • 1
    It seems that parts of your answer are missing, since neither code nor text allows anyone to calculate sunrise/sunset times. Dec 1, 2018 at 8:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.