[postgis-tickets] r15884 - ST_Angle function from R?\195?\169mi Cura
Regina Obe
lr at pcorp.us
Wed Oct 4 09:19:30 PDT 2017
Author: robe
Date: 2017-10-04 09:19:29 -0700 (Wed, 04 Oct 2017)
New Revision: 15884
Modified:
trunk/NEWS
trunk/doc/reference_measure.xml
trunk/postgis/lwgeom_functions_basic.c
trunk/postgis/postgis.sql.in
trunk/regress/lwgeom_regress.sql
trunk/regress/lwgeom_regress_expected
Log:
ST_Angle function from R?\195?\169mi Cura
Closes https://github.com/postgis/postgis/pull/97
Closes #3876
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2017-10-04 16:17:03 UTC (rev 15883)
+++ trunk/NEWS 2017-10-04 16:19:29 UTC (rev 15884)
@@ -1,8 +1,8 @@
PostGIS 2.5.0
2018/xx/xx
* New Features *
+ - #3876, ST_Angle function (Rémi Cura)
-
PostGIS 2.4.0
2017/09/30
Modified: trunk/doc/reference_measure.xml
===================================================================
--- trunk/doc/reference_measure.xml 2017-10-04 16:17:03 UTC (rev 15883)
+++ trunk/doc/reference_measure.xml 2017-10-04 16:19:29 UTC (rev 15884)
@@ -820,7 +820,73 @@
</refsection>
</refentry>
+
+ <refentry id="ST_Angle">
+ <refnamediv>
+ <refname>ST_Angle</refname>
+ <refpurpose>Returns the angle between 3 points, or between 2 vectors (4 points or 2 lines).</refpurpose>
+ </refnamediv>
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>float <function>ST_Angle</function></funcdef>
+ <paramdef><type>geometry </type><parameter>point1</parameter></paramdef>
+ <paramdef><type>geometry </type><parameter>point2</parameter></paramdef>
+ <paramdef><type>geometry </type><parameter>point3</parameter></paramdef>
+ <paramdef choice="opt"><type>geometry </type><parameter>point4</parameter></paramdef>
+ </funcprototype>
+ <funcprototype>
+ <funcdef>float <function>ST_Angle</function></funcdef>
+ <paramdef><type>geometry </type><parameter>line1</parameter></paramdef>
+ <paramdef><type>geometry </type><parameter>line2</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+ <refsection>
+ <title>Description</title>
+
+ <para> For 3 points, computes the angle measured clockwise of P1P2P3.
+ If input are 2 lines, get first and last point of the lines as 4 points.
+ For 4 points,compute the angle measured clockwise of P1P2,P3P4.
+ Results are always positive, between 0 and 2*Pi radians.
+
+ Uses azimuth of pairs or points.
+ </para>
+ <para>ST_Angle(P1,P2,P3) = ST_Angle(P2,P1,P2,P3)</para>
+ <para>Result is in radian and can be converted to degrees using a built-in PostgreSQL function degrees(), as shown in the example.</para>
+ <para>Availability: 2.5.0</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+ <para>Geometry Azimuth in degrees </para>
+<programlisting>
+ WITH rand AS (
+ SELECT s, random() * 2 * PI() AS rad1
+ , random() * 2 * PI() AS rad2
+ FROM generate_series(1,2,2) AS s
+ )
+ , points AS (
+ SELECT s, rad1,rad2, ST_MakePoint(cos1+s,sin1+s) as p1, ST_MakePoint(s,s) AS p2, ST_MakePoint(cos2+s,sin2+s) as p3
+ FROM rand
+ ,cos(rad1) cos1, sin(rad1) sin1
+ ,cos(rad2) cos2, sin(rad2) sin2
+ )
+ SELECT s, ST_AsText(ST_SnapToGrid(ST_MakeLine(ARRAY[p1,p2,p3]),0.001)) AS line
+ , degrees(ST_Angle(p1,p2,p3)) as computed_angle
+ , round(degrees(2*PI()-rad2 -2*PI()+rad1+2*PI()))::int%360 AS reference
+ , round(degrees(2*PI()-rad2 -2*PI()+rad1+2*PI()))::int%360 AS reference
+ FROM points ;
+
+1 | line | computed_angle | reference
+------------------+------------------
+1 | LINESTRING(1.511 1.86,1 1,0.896 0.005) | 155.27033848688 | 155
+
+</programlisting>
+ </refsection>
+ </refentry>
+
<refentry id="ST_Centroid">
<refnamediv>
<refname>ST_Centroid</refname>
Modified: trunk/postgis/lwgeom_functions_basic.c
===================================================================
--- trunk/postgis/lwgeom_functions_basic.c 2017-10-04 16:17:03 UTC (rev 15883)
+++ trunk/postgis/lwgeom_functions_basic.c 2017-10-04 16:19:29 UTC (rev 15884)
@@ -102,6 +102,7 @@
Datum LWGEOM_asEWKT(PG_FUNCTION_ARGS);
Datum LWGEOM_hasBBOX(PG_FUNCTION_ARGS);
Datum LWGEOM_azimuth(PG_FUNCTION_ARGS);
+Datum LWGEOM_angle(PG_FUNCTION_ARGS);
Datum LWGEOM_affine(PG_FUNCTION_ARGS);
Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS);
Datum optimistic_overlap(PG_FUNCTION_ARGS);
@@ -2420,6 +2421,131 @@
PG_RETURN_FLOAT8(result);
}
+/**
+ * Compute the angle defined by 3 points or the angle between 2 vectors
+ * defined by 4 points
+ * given Point geometries.
+ * @return NULL on exception (same point).
+ * Return radians otherwise (always positive).
+ */
+PG_FUNCTION_INFO_V1(LWGEOM_angle);
+Datum LWGEOM_angle(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED * seri_geoms[4];
+ LWGEOM *geom_unser;
+ LWPOINT *lwpoint;
+ POINT2D points[4];
+ double az1,az2 ;
+ double result;
+ int srids[4];
+ int i = 0 ;
+ int j = 0;
+ int err_code = 0;
+ int n_args = PG_NARGS();
+
+ /* no deserialize, checking for common error first*/
+ for(i=0; i<n_args; i++)
+ {
+ seri_geoms[i] = PG_GETARG_GSERIALIZED_P(i);
+ if (gserialized_is_empty(seri_geoms[i]) )
+ {/* empty geom */
+ if (i==3)
+ {
+ n_args = 3 ;
+ }
+ else
+ {
+ err_code = 1 ;
+ break ;
+ }
+ } else
+ {
+ if(gserialized_get_type(seri_geoms[i]) != POINTTYPE)
+ {/* geom type */
+ err_code = 2 ;
+ break;
+ }
+ else
+ {
+ srids[i] = gserialized_get_srid(seri_geoms[i]) ;
+ if(srids[0] != srids[i])
+ {/* error on srid*/
+ err_code = 3 ;
+ break;
+ }
+ }
+ }
+ }
+ if (err_code >0)
+ switch (err_code){
+ default: /*always executed*/
+ for (j=0;j<=i;j++)
+ PG_FREE_IF_COPY(seri_geoms[j], j);
+
+ case 1:
+ lwpgerror("Empty geometry");
+ PG_RETURN_NULL() ;
+ break;
+
+ case 2:
+ lwpgerror("Argument must be POINT geometries");
+ PG_RETURN_NULL();
+ break;
+
+ case 3:
+ lwpgerror("Operation on mixed SRID geometries");
+ PG_RETURN_NULL();
+ break;
+ }
+ /* extract points */
+ for(i=0; i<n_args; i++)
+ {
+ geom_unser = lwgeom_from_gserialized(seri_geoms[i]) ;
+ lwpoint = lwgeom_as_lwpoint(geom_unser);
+ if (!lwpoint)
+ {
+ for (j=0;j<n_args;j++)
+ PG_FREE_IF_COPY(seri_geoms[j], j);
+ lwpgerror("Error unserializing geometry");
+ PG_RETURN_NULL() ;
+ }
+
+ if ( ! getPoint2d_p(lwpoint->point, 0, &points[i]) )
+ {
+ /* // can't free serialized geom, it might be needed by lw
+ for (j=0;j<n_args;j++)
+ PG_FREE_IF_COPY(seri_geoms[j], j); */
+ lwpgerror("Error extracting point");
+ PG_RETURN_NULL();
+ }
+ /* lwfree(geom_unser);don't do, lw may rely on this memory
+ lwpoint_free(lwpoint); dont do , this memory is needed ! */
+ }
+ /* // can't free serialized geom, it might be needed by lw
+ for (j=0;j<n_args;j++)
+ PG_FREE_IF_COPY(seri_geoms[j], j); */
+
+ /* compute azimuth for the 2 pairs of points
+ * note that angle is not defined identically for 3 points or 4 points*/
+ if (n_args == 3)
+ {/* we rely on azimuth to complain if points are identical */
+ if ( ! azimuth_pt_pt(&points[0], &points[1], &az1) )
+ PG_RETURN_NULL();
+ if ( ! azimuth_pt_pt(&points[2], &points[1], &az2) )
+ PG_RETURN_NULL();
+ } else
+ {
+ if ( ! azimuth_pt_pt(&points[0], &points[1], &az1) )
+ PG_RETURN_NULL();
+ if ( ! azimuth_pt_pt(&points[2], &points[3], &az2) )
+ PG_RETURN_NULL();
+ }
+ result = az2-az1 ;
+ result += (result<0) * 2 * M_PI ; /* we dont want negative angle*/
+ PG_RETURN_FLOAT8(result);
+}
+
+
/*
* optimistic_overlap(Polygon P1, Multipolygon MP2, double dist)
* returns true if P1 overlaps MP2
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2017-10-04 16:17:03 UTC (rev 15883)
+++ trunk/postgis/postgis.sql.in 2017-10-04 16:19:29 UTC (rev 15884)
@@ -1290,7 +1290,15 @@
RETURNS float8
AS 'MODULE_PATHNAME', 'LWGEOM_azimuth'
LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL;
+
+-- Availability: 2.5.0
+CREATE OR REPLACE FUNCTION ST_Angle(pt1 geometry, pt2 geometry, pt3 geometry, pt4 geometry default 'POINT EMPTY'::geometry)
+ RETURNS float8
+ AS 'MODULE_PATHNAME', 'LWGEOM_angle'
+ LANGUAGE 'c' IMMUTABLE STRICT;
+
+
------------------------------------------------------------------------
-- MISC
------------------------------------------------------------------------
@@ -5974,6 +5982,17 @@
AS $$SELECT @extschema at ._ST_AsX3D(3,$1,$2,$3,'');$$
LANGUAGE 'sql' IMMUTABLE _PARALLEL;
+
+-----------------------------------------------------------------------
+-- ST_Angle
+-----------------------------------------------------------------------
+-- Availability: 2.3.0
+-- has to be here because need ST_StartPoint
+CREATE OR REPLACE FUNCTION ST_Angle(line1 geometry, line2 geometry)
+ RETURNS float8 AS 'SELECT ST_Angle(St_StartPoint($1), ST_EndPoint($1), St_StartPoint($2), ST_EndPoint($2))'
+ LANGUAGE 'sql' IMMUTABLE STRICT;
+
+
-- make views and spatial_ref_sys public viewable --
GRANT SELECT ON TABLE geography_columns TO public;
GRANT SELECT ON TABLE geometry_columns TO public;
Modified: trunk/regress/lwgeom_regress.sql
===================================================================
--- trunk/regress/lwgeom_regress.sql 2017-10-04 16:17:03 UTC (rev 15883)
+++ trunk/regress/lwgeom_regress.sql 2017-10-04 16:19:29 UTC (rev 15884)
@@ -135,3 +135,37 @@
SELECT 'BoundingDiagonal6', ST_AsEwkt(ST_BoundingDiagonal(
'SRID=3857;POLYGON M EMPTY'::geometry
));
+
+
+--- ST_Azimuth
+SELECT 'ST_Azimuth_regular' , round(ST_Azimuth(geom1,geom2)::numeric,4)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, CAST('POINT(1 0)' AS geometry) AS geom2 ;
+SELECT 'ST_Azimuth_same_point' , ST_Azimuth(geom1,geom1)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1 ;
+SELECT 'ST_Azimuth_mixed_srid' , ST_Azimuth(geom1,geom2)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, ST_GeomFromText('POINT(1 0)',4326) AS geom2;
+SELECT 'ST_Azimuth_not_point' , ST_Azimuth(geom1,geom2)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, ST_GeomFromText('LINESTRING(1 0 ,2 0)',4326) AS geom2;
+SELECT 'ST_Azimuth_null_geom' , ST_Azimuth(geom1,geom2)
+FROM CAST('POINT(0 1)' AS geometry) AS geom1, ST_GeomFromText('EMPTY') AS geom2;
+
+--- ST_Angle( points)
+SELECT 'ST_Angle_4_pts', St_Angle(p1,p2,p3,p4)
+ FROM ST_GeomFromtext('POINT(0 1)') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+ , ST_GeomFromtext('POINT(1 0)') AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+SELECT 'ST_Angle_4_pts', St_Angle(p1,p2,p3,p4)
+ FROM ST_GeomFromtext('POINT(2 0)') AS p1, ST_GeomFromtext('POINT(1 0)') AS p2
+ , ST_GeomFromtext('POINT(1 -1)') AS p3, ST_GeomFromtext('POINT(0 0)') AS p4 ;
+SELECT 'ST_Angle_3_pts', St_Angle(p1,p2,p3)
+ FROM ST_GeomFromtext('POINT(0 1)') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+ , ST_GeomFromtext('POINT(1 0)') AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+SELECT 'ST_Angle_mixed_srid', St_Angle(p1,p2,p3,p4)
+ FROM ST_GeomFromtext('POINT(0 1)') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+ , ST_GeomFromtext('POINT(1 0)',4326) AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+SELECT 'ST_Angle_empty' , St_Angle(p1,p2,p3,p4)
+ FROM ST_GeomFromtext('POINT EMPTY') AS p1, ST_GeomFromtext('POINT(0 0)') AS p2
+ , ST_GeomFromtext('POINT(1 0)',4326) AS p3, ST_GeomFromtext('POINT(2 0)') AS p4 ;
+--- ST_Angle( lines)
+SELECT 'ST_Angle_2_lines', St_Angle(l1,l2)
+ FROM ST_GeomFromtext('LINESTRING(0 1,0 0)') AS l1
+ , ST_GeomFromtext('LINESTRING(1 0, 2 0)') AS l2 ;
Modified: trunk/regress/lwgeom_regress_expected
===================================================================
--- trunk/regress/lwgeom_regress_expected 2017-10-04 16:17:03 UTC (rev 15883)
+++ trunk/regress/lwgeom_regress_expected 2017-10-04 16:19:29 UTC (rev 15884)
@@ -21,3 +21,14 @@
BoundingDiagonal4|SRID=3857;LINESTRING(-1 -2 -8 2,1 2 3 9)
BoundingDiagonal5|SRID=3857;LINESTRINGM(4 4 0,5 4 1)
BoundingDiagonal6|SRID=3857;LINESTRINGM EMPTY
+ST_Azimuth_regular|2.3562
+ST_Azimuth_same_point|
+ERROR: Operation on mixed SRID geometries
+ERROR: Argument must be POINT geometries
+ERROR: parse error - invalid geometry
+ST_Angle_4_pts|4.71238898038469
+ST_Angle_4_pts|0.785398163397448
+ST_Angle_3_pts|1.5707963267949
+ERROR: Operation on mixed SRID geometries
+ERROR: Empty geometry
+ST_Angle_2_lines|4.71238898038469
More information about the postgis-tickets
mailing list