Reputation: 1861
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.
Upvotes: 24
Views: 44239
Reputation: 496
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
Upvotes: 19
Reputation: 1626
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
;
Upvotes: 6
Reputation: 1
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
Upvotes: -3
Reputation: 117380
Take a look at these links:
Upvotes: 29
Reputation: 1180
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.
Upvotes: 12
Reputation: 37
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]
Upvotes: 2