[postgis-users] Solar azimuth calculator
Imre Samu
pella.samu at gmail.com
Sun Apr 18 06:33:02 PDT 2021
> Can this be added as a native function to Postgresql or postGIS? I'm
getting a lot of requests for this data.
IMHO: There are multiple other alternatives:
via PL/Python3
-
https://pvlib-python.readthedocs.io/en/stable/_modules/pvlib/solarposition.html
- https://pysolar.readthedocs.io/en/latest/
via PL/R https://rdrr.io/search?q=solar%20azimuth
via plv8 ( V8 Engine Javascript Procedural Language add-on for PostgreSQL )
- https://www.npmjs.com/package/suncalc
- ...
via Rust ( https://github.com/zombodb/pgx Build Postgres Extensions with
Rust! )
- https://github.com/frehberg/spa-rs ( Solar Position Algorithm
implementation in Rust )
- ...
Imre
Alexander Gataric <gataric at usa.net> ezt írta (időpont: 2021. ápr. 18., V,
15:01):
> Can this be added as a native function to Postgresql or postGIS? I'm
> getting a lot of requests for this data.
>
> Thanks
> Alex
>
> Get BlueMail for Android <http://www.bluemail.me/r?b=16696>
> On Apr 17, 2021, at 1:57 PM, Alexander Gataric <gataric at usa.net> wrote:
>>
>> Thanks!
>>
>> Get BlueMail for Android <http://www.bluemail.me/r?b=16696>
>> On Apr 17, 2021, at 9:12 AM, Bruce Rindahl < bruce.rindahl at gmail.com>
>> wrote:
>>>
>>> Here it is in psql. Doesn't need PostGIS. Fun little exercise.
>>>
>>>
>>> CREATE OR REPLACE FUNCTION public.solarazelev(
>>>
>>> obs_time timestamp with time zone,
>>>
>>> latitude numeric,
>>>
>>> longitude numeric,
>>>
>>> alt numeric)
>>>
>>> RETURNS TABLE(azimuth numeric, elevation numeric)
>>>
>>> LANGUAGE 'plpgsql'
>>>
>>> AS $BODY$
>>>
>>> DECLARE copyright varchar :=
>>>
>>> ' Programed by Darin C.Koblick 2 / 17 / 2009
>>>
>>> Darin C.Koblick 4 / 16 / 2013 Vectorized for Speed
>>>
>>> Allow for MATLAB Datevec input
>>> in
>>>
>>> addition to a UTC string.
>>>
>>> Cleaned up comments and code to
>>>
>>> avoid warnings in MATLAB editor.
>>>
>>> Kevin Godden 9/1/2020 Ported from Matlab to C++,
>>> tried to change as little as possible.
>>>
>>> this is a non-vectorised port.
>>>
>>> Bruce Rindahl 4/17/2021 Ported to PostgreSQL
>>>
>>>
>>> --------------------------------------------------------------------------
>>>
>>> Copyright(c) 2010, Darin Koblick
>>>
>>> All rights reserved.
>>>
>>>
>>>
>>> Redistribution and use in source and binary forms, with or without
>>>
>>> modification, are permitted provided that the following conditions are
>>> met :
>>>
>>>
>>>
>>> *Redistributions of source code must retain the above copyright notice,
>>> this
>>>
>>> list of conditions and the following disclaimer.
>>>
>>>
>>>
>>> * Redistributions in binary form must reproduce the above copyright
>>> notice,
>>>
>>> this list of conditions and the following disclaimer in the documentation
>>>
>>> and/or other materials provided with the distribution
>>>
>>>
>>>
>>> THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
>>> IS"
>>>
>>> AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
>>>
>>> IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
>>> PURPOSE ARE
>>>
>>> DISCLAIMED.IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
>>> LIABLE
>>>
>>> FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
>>> CONSEQUENTIAL
>>>
>>> DAMAGES(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
>>>
>>> SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
>>> HOWEVER
>>>
>>> CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
>>> LIABILITY,
>>>
>>> OR TORT(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
>>> USE
>>>
>>> OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.';
>>>
>>>
>>>
>>> DECLARE M_PI numeric := pi();
>>>
>>>
>>>
>>> DECLARE jd numeric := to_char(obs_time,'J')::numeric;
>>>
>>>
>>>
>>> DECLARE d numeric := jd - 2451543.5;
>>>
>>>
>>>
>>> -- Keplerian Elements for the Sun(geocentric)
>>>
>>> DECLARE w numeric := 282.9404 + 4.70935e-5*d; -- (longitude of
>>> perihelion degrees)
>>>
>>> -- a = 1.000000; % (mean distance, a.u.)
>>>
>>> DECLARE e numeric := 0.016709 - 1.151e-9*d; -- (eccentricity)
>>>
>>> DECLARE M numeric := mod(356.0470 + 0.9856002585*d, 360.0); --
>>> (mean anomaly degrees)
>>>
>>> DECLARE L numeric := w + M; -- (Sun's mean longitude degrees)
>>>
>>>
>>>
>>> DECLARE oblecl numeric := 23.4393 - 3.563e-7*d; -- (Sun's obliquity
>>> of the ecliptic)
>>>
>>> -- auxiliary angle
>>>
>>> DECLARE EE numeric := M + (180 / M_PI)*e*sin(M*(M_PI / 180))*(1 +
>>> e*cos(M*(M_PI / 180)));
>>>
>>>
>>>
>>> -- rectangular coordinates in the plane of the ecliptic(x axis
>>> toward perhilion)
>>>
>>> DECLARE x numeric := cos(EE*(M_PI / 180)) - e;
>>>
>>> DECLARE y numeric := sin(EE*(M_PI / 180))*sqrt(1 - pow(e, 2));
>>>
>>>
>>>
>>> -- find the distance and true anomaly
>>>
>>> DECLARE r numeric := sqrt(pow(x,2) + pow(y,2));
>>>
>>> -- Is this y,x or x,y? opposite from Excel - use C++ version
>>>
>>> DECLARE v numeric := atan2(y, x)*(180 / M_PI);
>>>
>>>
>>>
>>> -- find the longitude of the sun
>>>
>>> DECLARE lon numeric := v + w;
>>>
>>>
>>>
>>> -- compute the ecliptic rectangular coordinates
>>>
>>> DECLARE xeclip numeric := r*cos(lon*(M_PI / 180));
>>>
>>> DECLARE yeclip numeric := r*sin(lon*(M_PI / 180));
>>>
>>> DECLARE zeclip numeric := 0.0;
>>>
>>> --rotate these coordinates to equitorial rectangular coordinates
>>>
>>> DECLARE xequat numeric := xeclip;
>>>
>>>
>>>
>>> DECLARE yequat numeric := yeclip*cos(oblecl*(M_PI / 180)) + zeclip
>>> * sin(oblecl*(M_PI / 180));
>>>
>>>
>>>
>>> DECLARE zequat numeric := yeclip*sin(23.4406*(M_PI / 180)) + zeclip
>>> * cos(oblecl*(M_PI / 180));
>>>
>>> -- convert equatorial rectangular coordinates to RA and Decl:
>>>
>>> DECLARE RR numeric := sqrt(pow(xequat, 2) + pow(yequat, 2) +
>>> pow(zequat, 2)) - (Alt / 149598000); -- roll up the altitude correction
>>>
>>> DECLARE RA numeric := atan2(yequat, xequat)*(180 / M_PI);
>>>
>>>
>>>
>>> DECLARE delta numeric := asin(zequat / RR)*(180 / M_PI);
>>>
>>>
>>>
>>> DECLARE UTH numeric := date_part('hour',obs_time)
>>>
>>> + date_part('minute',obs_time)/60
>>>
>>> + date_part('second',obs_time) /3600
>>>
>>> - extract( timezone from now())/3600;
>>>
>>> -- Calculate local siderial time
>>>
>>> DECLARE GMST0 numeric := mod(L + 180, 360.0) / 15;
>>>
>>>
>>>
>>> DECLARE SIDTIME numeric := GMST0 + UTH + Longitude / 15;
>>>
>>>
>>>
>>> -- Replace RA with hour angle HA
>>>
>>> DECLARE HA numeric := (SIDTIME*15 - RA);
>>>
>>>
>>>
>>> -- convert to rectangular coordinate system
>>>
>>> DECLARE XX numeric := cos(HA*(M_PI / 180))*cos(delta*(M_PI / 180));
>>>
>>>
>>>
>>> DECLARE YY numeric := sin(HA*(M_PI / 180))*cos(delta*(M_PI / 180));
>>>
>>> DECLARE z numeric := sin(delta*(M_PI / 180));
>>>
>>>
>>>
>>> -- rotate this along an axis going east - west.
>>>
>>> DECLARE xhor numeric := XX*cos((90 - Latitude)*(M_PI / 180)) -
>>> z*sin((90 - Latitude)*(M_PI / 180));
>>>
>>>
>>>
>>> DECLARE yhor numeric := YY;
>>>
>>> DECLARE zhor numeric := XX*sin((90 - Latitude)*(M_PI / 180)) +
>>> z*cos((90 - Latitude)*(M_PI / 180));
>>>
>>>
>>>
>>> -- Find the h and AZ
>>>
>>> DECLARE Az numeric := round((atan2(yhor, xhor)*(180 / M_PI) +
>>> 180)::numeric,3);
>>>
>>> DECLARE El numeric := round((asin(zhor)*(180 / M_PI))::numeric,3);
>>>
>>> BEGIN
>>>
>>> RETURN QUERY select Az,El;
>>>
>>> END
>>>
>>> $BODY$;
>>>
>>>
>>>
>>>
>>> On Thu, Apr 15, 2021, 1:46 PM Bruce Rindahl < bruce.rindahl at gmail.com>
>>> wrote:
>>>
>>>> Here is the code for it in c:
>>>>
>>>>
>>>> https://www.ridgesolutions.ie/index.php/2020/01/14/c-code-to-estimate-solar-azimuth-and-elevation-given-gps-position-and-time/
>>>>
>>>>
>>> ------------------------------
>>>
>>> postgis-users mailing list
>>> postgis-users at lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/postgis-users
>>>
>>> ------------------------------
>>
>> postgis-users mailing list
>> postgis-users at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/postgis-users
>>
>> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210418/51387c6e/attachment.html>
More information about the postgis-users
mailing list