在PostgreSQL中计算给定GPS坐标的日出和日落时间。

Question

在PostgreSQL中计算给定GPS坐标的日出和日落时间。

pythonpostgresqlgps

24

我想将PostgreSQL表中的timestamp数据类型分类，以确定它们是否可视为“白天”或“黑夜”。换句话说，我希望能够精确地计算出在特定GPS位置下的日出和日落时间。

我知道plpgsql和plpython。

- zaadeh

1

根据被接受的答案，似乎可以将标题从“在PostgreSQL内部”更改为“在Python中”。这个解释正确吗？ - reducing activity

6个回答

19

我知道这个很老了，但由于我找不到快速的解决方案，所以我想分享一下。这里使用了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']

至少在我所在的地方，它似乎准确到几分钟。为了更高的精度，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

- oortCloud

太棒了！非常有用的小类。 - Nico Coallier

链接已失效。不过，它可以在 Wayback Machine（https://web.archive.org/web/20161202180207/http://williams.best.vwh.net:80/sunrise_sunset_algorithm.htm）下进行归档。 - Muhammad Yasirroni

当前原作者的链接为https://edwilliams.org/sunrise_sunset_algorithm.htm。 - Muhammad Yasirroni

11

使用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)

给出：

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

你可以使用这个作为起点，编写你自己的postgres（或postgis）函数，但要使用 plpython 而不是 plr。

- jonas

对于任何位置，您都可以调用 Astral().sun_utc(date.today(), latitude, longitude) 或其他变量来获取当地时间、不同的黄昏时段等信息。 - rcoup

1

自1.6版本以来，astral API已经发生了很大的变化，但仍然是一个不错的工具-请查看上面链接的文档以获取示例用法。 - FObersteiner

6

我知道这是一个古老的问题，但我一直需要在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;

使用示例：

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

- Adam Charnock

2

我使用这个工具来计算日出、日落、黎明和黄昏时间。

只需将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]

- Chris

5

请至少包含你的导入项。 - Daniel Moodie

-3

我正在使用Sun.py。今天我得到了minutes = 60，UT = 12.9979740551的值。

类Sun：

#10. Return

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

#add this after min calculate

        if min == 60:
            hr += 1
            min = 0

- Aleksey Konovalov

1

似乎你的回答缺少部分内容，因为无论是代码还是文本都不能让任何人计算日出/日落时间。 - reducing activity

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Roman Pekar · Accepted Answer

看看这些链接：

使用Python计算日出和日落时间;
Skyfield 项目 (PyEphem的新版本)
PyEphem 项目;
astral 项目;