Python 计算 PostgreSQL 中给定 GPS 坐标的日出和日落时间
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/19615350/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Calculate sunrise and sunset times for a given GPS coordinate within PostgreSQL
提问by zaadeh
I want to classify timestampdata 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.
我想对timestampPostgreSQL 表中的数据类型进行分类，看它们是“白天”还是“晚上”。换句话说，给定特定的 GPS 位置，我希望能够准确计算日出和日落时间。
I know plpgsql and plpython.
我知道 plpgsql 和 plpython。
采纳答案by Roman Pekar
Take a look at these links:
看看这些链接：
- Calulating sunrise and sunset in Python;
- Skyfieldproject (new incarnation of PyEphem)
- PyEphemproject;
- astralproject;
- 用 Python 计算日出和日落；
- Skyfield项目（PyEphem 的新化身）
- PyEphem项目；
- 星体项目；
回答by Chris
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.
只需要用您的坐标替换 X。请注意，返回的时间是 UTC，因此您需要根据需要添加您的特定时区小时数。
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]
回答by oortCloud
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.
我知道这已经过时了，但我想我会分享，因为我找不到快速解决方案。这使用了 Sun 类（见下文），我按照此链接构建了该类。
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.
它似乎准确到几分钟内，至少我住的地方。为了获得更高的精度，calcSunTime() 方法中的 zenith 参数可以使用微调。有关更多信息，请参阅上面的链接。
# 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
回答by jonas
Use Astral(current version 1.6). The first example in the documentationshows the calculation of sunrise and sunset for a given location. A simpler example with custom latitude and longitude would be:
使用Astral（当前版本 1.6）。该文档中的第一实施例示出了日出和日落的对于给定的位置的计算。自定义纬度和经度的一个更简单的例子是：
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 thisas a starting point for writing your own postgres (or postgis) functions using plpython instead of plr.
然后，您可以将其用作使用plpython 而不是 plr 编写自己的 postgres（或 postgis）函数的起点。
回答by Aleksey Konovalov
I'm using Sun.py. Today I got the value of minutes = 60, UT = 12.9979740551
我正在使用 Sun.py。今天，我得到的价值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
回答by Adam Charnock
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:
所以我知道这是一个古老的问题，但我一直需要在 Postgres 中实际计算它。因此，我将 oortCloud 的答案移植到了 PL/pgSQL。也许它对某人有用：
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
;
