[postgis-users] Reduce precision using ST_union with gridSize

[Palacios García Cuevas, Elena]

Hello, 

I am trying to simplify a polygon by reducing the precision using ST_Union with the gridSize parameter ( I am using  postgis 3.1  and GEOS 3.9.1). I have read that it is possible to do that by a union of a geometry with an empty geometry (http://blog.cleverelephant.ca/2020/12/waiting-postgis-31-3.html).   However, does not works for me ( maybe is my fault because I am new with postgis). I have tried as follows:

CREATE TABLE public.error_union_empty

(

    geom geometry(MultiPolygon),

    gid integer NOT NULL DEFAULT nextval('fallo_gid_seq'::regclass),

    CONSTRAINT error_union_empty_pkey PRIMARY KEY (gid)

);

INSERT INTO error_union_empty

(SELECT ST_union(geom,ST_GeomFromText('POLYGON EMPTY'), 0.001152) geom

FROM f_polygon);

But the result is an empty table.

Then I tried to do the same but by doing the union with the polygon himself, and this approach seems to work ok.

CREATE TABLE public.error_union_empty2

(

    geom geometry(MultiPolygon),

    gid integer NOT NULL DEFAULT nextval('fallo_gid_seq'::regclass),

    CONSTRAINT error_union_empty2_pkey PRIMARY KEY (gid)

);

insert into error_union_empty2

select  ST_CollectionExtract(ST_union(f_polygon.geom,f_polygon.geom, 0.001152)) geom

from f_polygon;

Why is not working the first approach? Have I missed something in the query?

Any help or comments are welcome.

Does anyone knows if  is expected to have a simplify with gridsize in the future?

Thanks!

---

Este correo electrónico y, en su caso, cualquier fichero anexo al mismo, contiene información de carácter confidencial exclusivamente dirigida a su destinatario o destinatarios. Si no es vd. el destinatario indicado, queda notificado que la lectura, utilización, divulgación y/o copia sin autorización está prohibida en virtud de la legislación vigente. En el caso de haber recibido este correo electrónico por error, se ruega notificar inmediatamente esta circunstancia mediante reenvío a la dirección electrónica del remitente.
Evite imprimir este mensaje si no es estrictamente necesario.

This email and any file attached to it (when applicable) contain(s) confidential information that is exclusively addressed to its recipient(s). If you are not the indicated recipient, you are informed that reading, using, disseminating and/or copying it without authorisation is forbidden in accordance with the legislation in effect. If you have received this email by mistake, please immediately notify the sender of the situation by resending it to their email address.
Avoid printing this message if it is not absolutely necessary.

[Martin Davis]

The issue with ST_Union not working with an empty polygon certainly sounds like a bug, and we should look into that. If you can file a ticket that would be great.   (It would be helpful to have the polygon data which is causing the problem.  Is it valid?)

In any case now there is the function ST_ReducePrecision [1] which is specifically designed to do what you want.  (It uses most of the same machinery as ST_Union, but should avoid whatever the issue with ST_Union is).   

[1] https://postgis.net/docs/manual-3.1/ST_ReducePrecision.html

_______________________________________________
postgis-users mailing list
postgi...@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[Palacios García Cuevas, Elena]

Thank you Martin.  With ST_ReducePrecision I got the expected result. The polygon, although is a complex one, is valid. I will file a ticket with a sample for the ST_Union issue.

---

De: postgis-users <postgis-us...@lists.osgeo.org> en nombre de Martin Davis <mtnc...@gmail.com>
Enviado: martes, 6 de abril de 2021 20:46
Para: PostGIS Users Discussion <postgi...@lists.osgeo.org>
Asunto: Re: [postgis-users] Reduce precision using ST_union with gridSize

 

[Alexander Gataric]

Does anyone have the SQL for PostGIS to calculate the solar azimuth and elevation for a time and point?

Thanks

Get BlueMail for Android

[Bruce Rindahl]

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/

[Bruce Rindahl]

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$;

[Alexander Gataric]

Thanks!

Get BlueMail for Android

[Alexander Gataric]

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

[Imre Samu]

> 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

_______________________________________________

[Bruce Rindahl]

Why do you need it as a native function?  Even if it was it wouldn't be available until PG 15/Postgis 3.3 at the earliest and then you would force everyone to upgrade.  I think this function would work on PG 8.

Remember PG SQL is Turing complete.  This means anything you want to do can be done in pure SQL. The only benefit of a native function is speed/memory.

I think a better idea is to add it to the users wiki https://trac.osgeo.org/postgis/wiki/UsersWikiMain and promote that.

[Bruce Rindahl]

As an added bonus putting it in the wiki makes bug fixes easier.

The line:

 - extract( timezone from now())/3600;

should read:

 - extract( timezone from obs_time)/3600;

;)

[Alexander Gataric]

The reason is that I need special permission to use this at my company since it is not part of the standard package. Technically there is no issue.

Get BlueMail for Android