[postgis-tickets] [SCM] PostGIS branch master updated. 3.3.0rc2-1091-gaeb6f6aa6

git at osgeo.org git at osgeo.org
Fri Jun 30 15:07:55 PDT 2023


This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".

The branch, master has been updated
       via  aeb6f6aa634542af74b3cc9b9792d77ad7798aab (commit)
      from  494e4f7c64558fee8335662144cc95ce4debe234 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit aeb6f6aa634542af74b3cc9b9792d77ad7798aab
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Fri Jun 30 15:06:25 2023 -0700

    Add extra signatures for ST_Project,
    
    * ST_Project(geom, distance, azimuth)
    * ST_Project(geom, geom, distance)
    * ST_Project(geog, geog, distance)
    
    In partial fulfillment of #5267, but still
    probably need an explicit line extender.

diff --git a/NEWS b/NEWS
index 1ae3f983d..14228e179 100644
--- a/NEWS
+++ b/NEWS
@@ -23,6 +23,7 @@ PostGIS 3.4.0dev
   - GH721, New window-based ST_ClusterWithinWin and ST_ClusterIntersectingWin (Paul Ramsey)
   - #5397, [address_standardizer] debug_standardize_address function (Regina Obe)
   - ST_LargestEmptyCircle, exposes extra semantics on circle finding (Martin Davis)
+  - ST_Project signature for geometry, and two-point signature (Paul Ramsey)
 
 * Enhancements *
   - #5194, do not update system catalogs from postgis_extensions_upgrade (Sandro Santilli)
@@ -35,6 +36,7 @@ PostGIS 3.4.0dev
   - Repair spatial planner stats to use computed selectivity for contains/within queries (Paul Ramsey)
   - GH734, Additional metadata on Proj installation in postgis_proj_version (Paul Ramsey)
   - #5177, allow building tools without PostgreSQL server headers (Sandro Santilli)
+  - ST_Project signature for geometry, and two-point signature (Paul Ramsey)
 
 * Bug Fix *
 
diff --git a/doc/reference_measure.xml b/doc/reference_measure.xml
index a9091ece9..1e94574c4 100644
--- a/doc/reference_measure.xml
+++ b/doc/reference_measure.xml
@@ -1879,36 +1879,67 @@ FROM ST_GeogFromText('MULTIPOLYGON(((-71.1044543107478 42.340674480411,-71.10445
 		  <refsynopsisdiv>
 			<funcsynopsis>
 			  <funcprototype>
-				<funcdef>geography <function>ST_Project</function></funcdef>
-
-				<paramdef><type>geography </type>
+				<funcdef>geometry <function>ST_Project</function></funcdef>
+				<paramdef><type>geometry </type>
 				<parameter>g1</parameter></paramdef>
 				<paramdef><type>float </type>
 				<parameter>distance</parameter></paramdef>
 				<paramdef><type>float </type>
 				<parameter>azimuth</parameter></paramdef>
 			  </funcprototype>
+
+              <funcprototype>
+                <funcdef>geometry <function>ST_Project</function></funcdef>
+                <paramdef><type>geometry </type>
+                <parameter>g1</parameter></paramdef>
+                <paramdef><type>geometry </type>
+                <parameter>g2</parameter></paramdef>
+                <paramdef><type>float </type>
+                <parameter>distance</parameter></paramdef>
+              </funcprototype>
+
+              <funcprototype>
+                <funcdef>geography <function>ST_Project</function></funcdef>
+                <paramdef><type>geography </type>
+                <parameter>g1</parameter></paramdef>
+                <paramdef><type>float </type>
+                <parameter>distance</parameter></paramdef>
+                <paramdef><type>float </type>
+                <parameter>azimuth</parameter></paramdef>
+              </funcprototype>
+
+              <funcprototype>
+                <funcdef>geography <function>ST_Project</function></funcdef>
+                <paramdef><type>geography </type>
+                <parameter>g1</parameter></paramdef>
+                <paramdef><type>geography </type>
+                <parameter>g2</parameter></paramdef>
+                <paramdef><type>float </type>
+                <parameter>distance</parameter></paramdef>
+              </funcprototype>
+
 			</funcsynopsis>
 		  </refsynopsisdiv>
 
 		  <refsection>
 			<title>Description</title>
 
-			<para>Returns a point projected from a start point along a geodesic using
-			a given distance and azimuth (bearing).
-			This is known as the direct geodesic problem.</para>
+			<para>Returns a point projected from a point along a geodesic using a given distance and azimuth (bearing). This is known as the direct geodesic problem.</para>
+            <para>The two-point version uses the path from the first to the second point to implicitly define the azimuth and uses the distance as before.</para>
 			<para>The distance is given in meters.  Negative values are supported.</para>
 			<para>The azimuth (also known as heading or bearing) is given in radians.
-			It is measured clockwise from true north (azimuth zero).
-			East is azimuth &#x03C0;/2 (90 degrees);
-			south is azimuth &#x03C0; (180 degrees);
-			west is azimuth 3&#x03C0;/2 (270 degrees).
-			Negative azimuth values and values greater than 2&#x03C0; (360 degrees) are supported.
-			</para>
-
+			It is measured clockwise from true north.</para>
+            <itemizedlist>
+                <listitem><para>North is azimuth zero (0 degrees)</para></listitem>
+                <listitem><para>East is azimuth &#x03C0;/2 (90 degrees)</para></listitem>
+                <listitem><para>South is azimuth &#x03C0; (180 degrees)</para></listitem>
+                <listitem><para>West is azimuth 3&#x03C0;/2 (270 degrees)</para></listitem>
+            </itemizedlist>
+            <para>Negative azimuth values and values greater than 2&#x03C0; (360 degrees) are supported.</para>
 
 			<para role="availability" conformance="2.0.0">Availability: 2.0.0</para>
 			<para role="enhanced" conformance="2.4.0">Enhanced: 2.4.0 Allow negative distance and non-normalized azimuth.</para>
+            <para role="enhanced" conformance="3.4.0">Enhanced: 3.4.0 Allow geometry arguments and two-point form omitting azimuth.</para>
 
 		  </refsection>
 
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index f79303a57..3bf5b4820 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1339,6 +1339,8 @@ extern LWPOINT* lwcompound_get_lwpoint(const LWCOMPOUND *lwcmp, uint32_t where);
 extern double ptarray_length_2d(const POINTARRAY *pts);
 extern int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring);
 extern int azimuth_pt_pt(const POINT2D *p1, const POINT2D *p2, double *ret);
+extern LWPOINT* lwpoint_project_lwpoint(const LWPOINT* lwpoint1, const LWPOINT* lwpoint2, double distance);
+extern LWPOINT* lwpoint_project(const LWPOINT* lwpoint1, double distance, double azimuth);
 extern int lwpoint_inside_circle(const LWPOINT *p, double cx, double cy, double rad);
 
 extern LWGEOM* lwgeom_reverse(const LWGEOM *lwgeom);
@@ -1744,6 +1746,11 @@ extern double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwge
 */
 extern LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth);
 
+/**
+* Calculate the location of a point on a spheroid, give a start point, end point and distance.
+*/
+extern LWPOINT* lwgeom_project_spheroid_lwpoint(const LWPOINT *from, const LWPOINT *to, const SPHEROID *spheroid, double distance);
+
 /**
 * Derive a new geometry with vertices added to ensure no vertex is more
 * than max_seg_length (in radians) from any other vertex.
diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index 4fa50ddc3..cc45fdeef 100644
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -253,6 +253,12 @@ int p4d_same(const POINT4D *p1, const POINT4D *p2);
 int p3d_same(const POINT3D *p1, const POINT3D *p2);
 int p2d_same(const POINT2D *p1, const POINT2D *p2);
 
+/*
+* Projections
+*/
+int project_pt(const POINT2D *P, double distance, double azimuth, POINT2D *R);
+int project_pt_pt(const POINT4D *A, const POINT4D *B, double distance, POINT4D *R);
+
 /*
 * Area calculations
 */
@@ -356,7 +362,9 @@ int ptarray_isccw(const POINTARRAY *pa);
 /*
 * Same
 */
+char ptarray_same2d(const POINTARRAY *pa1, const POINTARRAY *pa2);
 char ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2);
+char lwpoint_same2d(const LWPOINT *p1, const LWPOINT *p2);
 char lwpoint_same(const LWPOINT *p1, const LWPOINT *p2);
 char lwline_same(const LWLINE *p1, const LWLINE *p2);
 char lwpoly_same(const LWPOLY *p1, const LWPOLY *p2);
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 7d9f43f9b..89cf47589 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -2107,8 +2107,8 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
 	GEOGRAPHIC_POINT geo_source, geo_dest;
 	POINT4D pt_dest;
 	double x, y;
-	POINTARRAY *pa;
 	LWPOINT *lwp;
+	int has_z, has_m;
 
 	/* Normalize distance to be positive*/
 	if ( distance < 0.0 ) {
@@ -2129,6 +2129,8 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
 	/* Convert to ta geodetic point */
 	x = lwpoint_get_x(r);
 	y = lwpoint_get_y(r);
+	has_z = lwgeom_has_z(lwpoint_as_lwgeom(r));
+	has_m = lwgeom_has_m(lwpoint_as_lwgeom(r));
 	geographic_point_init(x, y, &geo_source);
 
 	/* Try the projection */
@@ -2140,20 +2142,26 @@ LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, dou
 	}
 
 	/* Build the output LWPOINT */
-	pa = ptarray_construct(0, 0, 1);
 	pt_dest.x = rad2deg(longitude_radians_normalize(geo_dest.lon));
 	pt_dest.y = rad2deg(latitude_radians_normalize(geo_dest.lat));
-	pt_dest.z = pt_dest.m = 0.0;
-	ptarray_set_point4d(pa, 0, &pt_dest);
-	lwp = lwpoint_construct(r->srid, NULL, pa);
+	pt_dest.z = has_z ? lwpoint_get_z(r) : 0.0;
+	pt_dest.m = has_m ? lwpoint_get_m(r) : 0.0;
+	lwp = lwpoint_make(r->srid, has_z, has_m, &pt_dest);
 	lwgeom_set_geodetic(lwpoint_as_lwgeom(lwp), LW_TRUE);
 	return lwp;
 }
 
+LWPOINT* lwgeom_project_spheroid_lwpoint(const LWPOINT *from, const LWPOINT *to, const SPHEROID *spheroid, double distance)
+{
+	double azimuth = lwgeom_azumith_spheroid(from, to, spheroid);
+	LWPOINT *lwp = lwgeom_project_spheroid(to, spheroid, distance, azimuth);
+	return lwp;
+}
+
 
 /**
 * Calculate a bearing (azimuth) given a source and destination point.
-* @param r - location of first point.
+https://accesd.desjardins.ca/coast* @param r - location of first point.
 * @param s - location of second point.
 * @param spheroid - spheroid definition.
 * @return azimuth - azimuth in radians.
diff --git a/liblwgeom/lwpoint.c b/liblwgeom/lwpoint.c
index 3f28c061e..d878c285f 100644
--- a/liblwgeom/lwpoint.c
+++ b/liblwgeom/lwpoint.c
@@ -266,6 +266,39 @@ lwpoint_same(const LWPOINT *p1, const LWPOINT *p2)
 	return ptarray_same(p1->point, p2->point);
 }
 
+/* check 2d coordinate equality  */
+char
+lwpoint_same2d(const LWPOINT *p1, const LWPOINT *p2)
+{
+	return ptarray_same2d(p1->point, p2->point);
+}
+
+LWPOINT *
+lwpoint_project_lwpoint(const LWPOINT* lwpoint1, const LWPOINT* lwpoint2, double distance)
+{
+	POINT4D p1, p2, p3;
+	int srid = lwgeom_get_srid((const LWGEOM*)lwpoint1);
+	int hasz = lwgeom_has_z((const LWGEOM*)lwpoint1);
+	int hasm = lwgeom_has_m((const LWGEOM*)lwpoint1);
+	lwpoint_getPoint4d_p(lwpoint1, &p1);
+	lwpoint_getPoint4d_p(lwpoint2, &p2);
+	project_pt_pt(&p1, &p2, distance, &p3);
+	return lwpoint_make(srid, hasz, hasm, &p3);
+}
+
+LWPOINT *
+lwpoint_project(const LWPOINT* lwpoint1, double distance, double azimuth)
+{
+	POINT4D p1, p2;
+	int srid = lwgeom_get_srid((const LWGEOM*)lwpoint1);
+	int hasz = lwgeom_has_z((const LWGEOM*)lwpoint1);
+	int hasm = lwgeom_has_m((const LWGEOM*)lwpoint1);
+	lwpoint_getPoint4d_p(lwpoint1, &p1);
+	lwpoint_getPoint4d_p(lwpoint1, &p2);
+	project_pt((POINT2D*)&p1, distance, azimuth, (POINT2D*)&p2);
+	return lwpoint_make(srid, hasz, hasm, &p2);
+}
+
 
 LWPOINT*
 lwpoint_force_dims(const LWPOINT *point, int hasz, int hasm, double zval, double mval)
diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index 6ee9e3426..d1e327cd2 100644
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -2466,3 +2466,49 @@ azimuth_pt_pt(const POINT2D *A, const POINT2D *B, double *d)
 	*d = fmod(2 * M_PI + M_PI / 2 - atan2(B->y - A->y, B->x - A->x), 2 * M_PI);
 	return LW_TRUE;
 }
+
+/**
+ * Azimuth is angle in radians from vertical axis
+ *
+ */
+int
+project_pt(const POINT2D *P, double distance, double azimuth, POINT2D *R)
+{
+	const double TWOPI = 2.0 * M_PI;
+	double slope;
+	/* Deal with azimuth out of (-360,360) range */
+	int orbits = floor(azimuth / TWOPI);
+	azimuth -= TWOPI * orbits;
+	/* Convert from azimuth to conventional slope */
+	slope = TWOPI - azimuth + M_PI_2;
+	if (slope > 0 && slope >  TWOPI) slope -= TWOPI;
+	if (slope < 0 && slope < -TWOPI) slope += TWOPI;
+
+	double dx = cos(slope) * distance;
+	double dy = sin(slope) * distance;
+	R->x += dx;
+	R->y += dy;
+	return LW_TRUE;
+}
+
+/**
+ * Azimuth is angle in radians from vertical axis
+ *
+ */
+int
+project_pt_pt(const POINT4D *A, const POINT4D *B, double distance, POINT4D *R)
+{
+	/* Convert from azimuth to conventional slope */
+	double len = distance2d_pt_pt((const POINT2D *)A, (const POINT2D *)B);
+	double prop = distance / len;
+	double dx = (B->x - A->x) * prop;
+	double dy = (B->y - A->y) * prop;
+	double dz = (B->z - A->z) * prop;
+	double dm = (B->m - A->m) * prop;
+	R->x += dx;
+	R->y += dy;
+	if (isfinite(dz)) R->z += dz;
+	if (isfinite(dm)) R->m += dm;
+	return LW_TRUE;
+}
+
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index d903d809d..6692485e7 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -497,6 +497,27 @@ ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2)
 	return LW_TRUE;
 }
 
+char
+ptarray_same2d(const POINTARRAY *pa1, const POINTARRAY *pa2)
+{
+	uint32_t i;
+
+	if ( FLAGS_GET_ZM(pa1->flags) != FLAGS_GET_ZM(pa2->flags) ) return LW_FALSE;
+	LWDEBUG(5,"dimensions are the same");
+
+	if ( pa1->npoints != pa2->npoints ) return LW_FALSE;
+	LWDEBUG(5,"npoints are the same");
+
+	for (i=0; i<pa1->npoints; i++)
+	{
+		if ( memcmp(getPoint_internal(pa1, i), getPoint_internal(pa2, i), sizeof(POINT2D)) )
+			return LW_FALSE;
+		LWDEBUGF(5,"point #%d is the same",i);
+	}
+
+	return LW_TRUE;
+}
+
 POINTARRAY *
 ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
 {
diff --git a/postgis/geography.sql.in b/postgis/geography.sql.in
index 6da02e208..9a61d7a7a 100644
--- a/postgis/geography.sql.in
+++ b/postgis/geography.sql.in
@@ -600,6 +600,13 @@ CREATE OR REPLACE FUNCTION ST_Project(geog geography, distance float8, azimuth f
 	LANGUAGE 'c' IMMUTABLE PARALLEL SAFE
 	_COST_MEDIUM;
 
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_Project(geog_from geography, geog_to geography, distance float8)
+	RETURNS geography
+	AS 'MODULE_PATHNAME','geography_project_geography'
+	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+	_COST_MEDIUM;
+
 -- Availability: 2.0.0
 CREATE OR REPLACE FUNCTION ST_Azimuth(geog1 geography, geog2 geography)
 	RETURNS float8
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 1911472c3..6bda676b3 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -1027,6 +1027,77 @@ Datum geography_project(PG_FUNCTION_ARGS)
 	PG_RETURN_POINTER(g_out);
 }
 
+/*
+** geography_project_geography(geog1, geog2, distance, azimuth)
+** returns point of projection given from/pt points,
+** azimuth in radians (bearing) and distance in meters
+*/
+PG_FUNCTION_INFO_V1(geography_project_geography);
+Datum geography_project_geography(PG_FUNCTION_ARGS)
+{
+	LWGEOM *lwgeom1, *lwgeom2;
+	LWPOINT *lwp1, *lwp2, *lwp3;
+	GSERIALIZED *g1, *g2, *g3;
+	double distance;
+	SPHEROID s;
+
+	/* Get our geometry object loaded into memory. */
+	g1 = PG_GETARG_GSERIALIZED_P(0);
+	g2 = PG_GETARG_GSERIALIZED_P(1);
+
+	if ( gserialized_get_type(g1) != POINTTYPE ||
+		 gserialized_get_type(g2) != POINTTYPE )
+	{
+		elog(ERROR, "ST_Project(geography) is only valid for point inputs");
+		PG_RETURN_NULL();
+	}
+
+	distance = PG_GETARG_FLOAT8(2); /* Distance in meters */
+
+	/* Handle the zero distance case */
+	if ( FP_EQUALS(distance, 0.0) )
+	{
+		PG_RETURN_POINTER(g2);
+	}
+
+	lwgeom1 = lwgeom_from_gserialized(g1);
+	lwgeom2 = lwgeom_from_gserialized(g2);
+
+	/* EMPTY things cannot be projected from */
+	if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
+	{
+		lwgeom_free(lwgeom1);
+		lwgeom_free(lwgeom2);
+		elog(ERROR, "ST_Project(geography) cannot project from an empty point");
+		PG_RETURN_NULL();
+	}
+
+	/* Initialize spheroid */
+	spheroid_init_from_srid(lwgeom_get_srid(lwgeom1), &s);
+
+	/* Calculate the length */
+	lwp1 = lwgeom_as_lwpoint(lwgeom1);
+	lwp2 = lwgeom_as_lwpoint(lwgeom2);
+	lwp3 = lwgeom_project_spheroid_lwpoint(lwp1, lwp2, &s, distance);
+
+	/* Something went wrong... */
+	if ( lwp3 == NULL )
+	{
+		elog(ERROR, "lwgeom_project_spheroid_lwpoint returned null");
+		PG_RETURN_NULL();
+	}
+
+	/* Clean up, but not all the way to the point arrays */
+	lwgeom_free(lwgeom1);
+	lwgeom_free(lwgeom2);
+	g3 = geography_serialize(lwpoint_as_lwgeom(lwp3));
+	lwpoint_free(lwp3);
+
+	PG_FREE_IF_COPY(g1, 0);
+	PG_FREE_IF_COPY(g2, 1);
+	PG_RETURN_POINTER(g3);
+}
+
 
 /*
 ** geography_azimuth(GSERIALIZED *g1, GSERIALIZED *g2)
diff --git a/postgis/lwgeom_functions_basic.c b/postgis/lwgeom_functions_basic.c
index 72e61bf77..4f22b76a9 100644
--- a/postgis/lwgeom_functions_basic.c
+++ b/postgis/lwgeom_functions_basic.c
@@ -105,6 +105,8 @@ Datum LWGEOM_setpoint_linestring(PG_FUNCTION_ARGS);
 Datum LWGEOM_asEWKT(PG_FUNCTION_ARGS);
 Datum LWGEOM_hasBBOX(PG_FUNCTION_ARGS);
 Datum LWGEOM_azimuth(PG_FUNCTION_ARGS);
+Datum geometry_project_direction(PG_FUNCTION_ARGS);
+Datum geometry_project_geometry(PG_FUNCTION_ARGS);
 Datum LWGEOM_angle(PG_FUNCTION_ARGS);
 Datum LWGEOM_affine(PG_FUNCTION_ARGS);
 Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS);
@@ -2573,6 +2575,82 @@ Datum LWGEOM_azimuth(PG_FUNCTION_ARGS)
 	PG_RETURN_FLOAT8(result);
 }
 
+
+/**
+ * Project a new point from a start point, direction and distance.
+ * ST_Project(geometry, distance, azimuth)
+ * Azimuth is measured in radians, clockwise from north.
+ * Distance is in SRID units.
+ * Geometry must be point.
+ */
+PG_FUNCTION_INFO_V1(geometry_project_direction);
+Datum geometry_project_direction(PG_FUNCTION_ARGS)
+{
+	GSERIALIZED *geom1, *geom2;
+	LWPOINT *lwpoint1, *lwpoint2;
+	LWGEOM *lwgeom1, *lwgeom2;
+	double distance, azimuth;
+
+	geom1 = PG_GETARG_GSERIALIZED_P(0);
+	distance = PG_GETARG_FLOAT8(1);
+	azimuth = PG_GETARG_FLOAT8(2);
+	lwgeom1 = lwgeom_from_gserialized(geom1);
+	lwpoint1 = lwgeom_as_lwpoint(lwgeom1);
+
+	if (!lwpoint1)
+		lwpgerror("Argument must be POINT geometry");
+
+	if (lwgeom_is_empty(lwgeom1))
+		PG_RETURN_NULL();
+
+	lwpoint2 = lwpoint_project(lwpoint1, distance, azimuth);
+	lwgeom2 = lwpoint_as_lwgeom(lwpoint2);
+	geom2 = geometry_serialize(lwgeom2);
+	PG_RETURN_POINTER(geom2);
+}
+
+
+/**
+ * Project a new point from a start point, direction and distance.
+ * ST_Project(geometry, distance, azimuth)
+ * Azimuth is measured in radians, clockwise from north.
+ * Distance is in SRID units.
+ * Geometry must be point.
+ */
+PG_FUNCTION_INFO_V1(geometry_project_geometry);
+Datum geometry_project_geometry(PG_FUNCTION_ARGS)
+{
+	GSERIALIZED *geom1, *geom2, *geom3;
+	LWPOINT *lwpoint1, *lwpoint2, *lwpoint3;
+	LWGEOM *lwgeom1, *lwgeom2, *lwgeom3;
+	double distance;
+
+	geom1 = PG_GETARG_GSERIALIZED_P(0);
+	geom2 = PG_GETARG_GSERIALIZED_P(1);
+	distance = PG_GETARG_FLOAT8(2);
+
+	lwgeom1 = lwgeom_from_gserialized(geom1);
+	lwpoint1 = lwgeom_as_lwpoint(lwgeom1);
+	lwgeom2 = lwgeom_from_gserialized(geom2);
+	lwpoint2 = lwgeom_as_lwpoint(lwgeom2);
+
+	if (!(lwpoint1 && lwpoint2))
+		lwpgerror("Arguments must be POINT geometries");
+
+	if (lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2))
+		PG_RETURN_NULL();
+
+	if (lwpoint_same2d(lwpoint1, lwpoint2))
+		PG_RETURN_POINTER(geom2);
+
+	lwpoint3 = lwpoint_project_lwpoint(lwpoint1, lwpoint2, distance);
+	lwgeom3 = lwpoint_as_lwgeom(lwpoint3);
+	geom3 = geometry_serialize(lwgeom3);
+
+	PG_RETURN_POINTER(geom3);
+}
+
+
 /**
  * Compute the angle defined by 3 points or the angle between 2 vectors
  * defined by 4 points
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 55d50df45..5e9275c74 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -1390,6 +1390,20 @@ CREATE OR REPLACE FUNCTION ST_azimuth(geom1 geometry, geom2 geometry)
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
 	_COST_LOW;
 
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_Project(geom1 geometry, distance float8, azimuth float8)
+	RETURNS geometry
+	AS 'MODULE_PATHNAME', 'geometry_project_direction'
+	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+	_COST_LOW;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_Project(geom1 geometry, geom2 geometry, distance float8)
+	RETURNS geometry
+	AS 'MODULE_PATHNAME', 'geometry_project_geometry'
+	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+	_COST_LOW;
+
 -- Availability: 2.5.0
 CREATE OR REPLACE FUNCTION ST_Angle(pt1 geometry, pt2 geometry, pt3 geometry, pt4 geometry default 'POINT EMPTY'::geometry)
 	RETURNS float8
diff --git a/regress/core/measures.sql b/regress/core/measures.sql
index 58072ca67..465578177 100644
--- a/regress/core/measures.sql
+++ b/regress/core/measures.sql
@@ -290,3 +290,19 @@ SELECT '#4328.3', ST_Intersects('TIN(((0 0,1 0,0 1,0 0)))'::geometry, 'TRIANGLE(
 SELECT '#4328.4', ST_Intersects('TIN(((0 0,1 0,0 1,0 0)))'::geometry, 'POLYGON((.1 .1, .2 .2, .2 .1, .1 .1))'::geometry), ST_3DIntersects('TIN(((0 0,1 0,0 1,0 0)))'::geometry, 'POLYGON((.1 .1, .2 .2, .2 .1, .1 .1))'::geometry);
 SELECT '#4328.5', ST_Intersects('TIN(((0 0,3 0,0 3,0 0)))'::geometry, 'CIRCULARSTRING(1.1 1.1, 1.2 1.2, 1.2 1.1)'::geometry), ST_3DIntersects('TIN(((0 0,3 0,0 3,0 0)))'::geometry, 'CIRCULARSTRING(1.1 1.1, 1.2 1.2, 1.2 1.1)'::geometry);
 SELECT '#4328.6', ST_Intersects('TIN(((0 0,3 0,0 3,0 0)))'::geometry, 'CURVEPOLYGON(CIRCULARSTRING(1.1 1.1, 1.2 1.2, 1.2 1.1, 1.2 1.2, 1.1 1.1))'::geometry), ST_3DIntersects('TIN(((0 0,3 0,0 3,0 0)))'::geometry, 'CURVEPOLYGON(CIRCULARSTRING(1.1 1.1, 1.2 1.2, 1.2 1.1, 1.2 1.2, 1.1 1.1))'::geometry);
+
+SELECT 'st_project.01', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 1, 0), 0.1), 2);
+SELECT 'st_project.02', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 1, pi()), 0.1), 2);
+SELECT 'st_project.03', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 1, pi()/2), 0.1), 2);
+SELECT 'st_project.04', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 1, 3*pi()/2), 0.1), 2);
+SELECT 'st_project.05', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 1, pi()/4), 0.001), 3);
+SELECT 'st_project.06', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 1, pi()+pi()/4), 0.001), 3);
+SELECT 'st_project.07', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 0, 0), 0.001), 3);
+
+SELECT 'st_project.11', ST_AsText(ST_SnapToGrid(ST_Project('POINT(1 0)'::geometry, 'POINT(0 0)'::geometry, 1), 0.1), 2);
+SELECT 'st_project.12', ST_AsText(ST_SnapToGrid(ST_Project('POINT(-1 0)'::geometry, 'POINT(0 0)'::geometry, 1), 0.1), 2);
+SELECT 'st_project.13', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 1)'::geometry, 'POINT(0 0)'::geometry, 1), 0.1), 2);
+SELECT 'st_project.14', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 -1)'::geometry, 'POINT(0 0)'::geometry, 1), 0.1), 2);
+SELECT 'st_project.15', ST_AsText(ST_SnapToGrid(ST_Project('POINT(1 1)'::geometry, 'POINT(0 0)'::geometry, 1), 0.001), 3);
+SELECT 'st_project.16', ST_AsText(ST_SnapToGrid(ST_Project('POINT(-1 -1)'::geometry, 'POINT(0 0)'::geometry, 1), 0.001), 3);
+SELECT 'st_project.17', ST_AsText(ST_SnapToGrid(ST_Project('POINT(0 0)'::geometry, 'POINT(0 0)'::geometry, 1), 0.001), 3);
diff --git a/regress/core/measures_expected b/regress/core/measures_expected
index 7ff9f1307..a9114589c 100644
--- a/regress/core/measures_expected
+++ b/regress/core/measures_expected
@@ -79,3 +79,17 @@ NOTICE:  One or both of the geometries is missing z-value. The unknown z-value w
 #4328.5|t|t
 NOTICE:  One or both of the geometries is missing z-value. The unknown z-value will be regarded as "any value"
 #4328.6|t|t
+st_project.01|POINT(0 1)
+st_project.02|POINT(0 -1)
+st_project.03|POINT(1 0)
+st_project.04|POINT(-1 0)
+st_project.05|POINT(0.707 0.707)
+st_project.06|POINT(-0.707 -0.707)
+st_project.07|POINT(0 0)
+st_project.11|POINT(-1 0)
+st_project.12|POINT(1 0)
+st_project.13|POINT(0 -1)
+st_project.14|POINT(0 1)
+st_project.15|POINT(-0.707 -0.707)
+st_project.16|POINT(0.707 0.707)
+st_project.17|POINT(0 0)

-----------------------------------------------------------------------

Summary of changes:
 NEWS                             |  2 ++
 doc/reference_measure.xml        | 57 ++++++++++++++++++++++-------
 liblwgeom/liblwgeom.h.in         |  7 ++++
 liblwgeom/liblwgeom_internal.h   |  8 +++++
 liblwgeom/lwgeodetic.c           | 20 +++++++----
 liblwgeom/lwpoint.c              | 33 +++++++++++++++++
 liblwgeom/measures.c             | 46 ++++++++++++++++++++++++
 liblwgeom/ptarray.c              | 21 +++++++++++
 postgis/geography.sql.in         |  7 ++++
 postgis/geography_measurement.c  | 71 ++++++++++++++++++++++++++++++++++++
 postgis/lwgeom_functions_basic.c | 78 ++++++++++++++++++++++++++++++++++++++++
 postgis/postgis.sql.in           | 14 ++++++++
 regress/core/measures.sql        | 16 +++++++++
 regress/core/measures_expected   | 14 ++++++++
 14 files changed, 375 insertions(+), 19 deletions(-)


hooks/post-receive
-- 
PostGIS


More information about the postgis-tickets mailing list