[postgis-users] Solar azimuth calculator
Alexander Gataric
gataric at usa.net
Sun Apr 18 06:00:41 PDT 2021
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
On Apr 17, 2021, 1:57 PM, at 1:57 PM, Alexander Gataric <gataric at usa.net> wrote:
>Thanks!
>
>Get BlueMail for Android
>
>On Apr 17, 2021, 9:12 AM, 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210418/cb081c48/attachment.html>
More information about the postgis-users
mailing list