self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-24 00:52:40 +00:00
Code Issues Packages Projects Releases Wiki Activity

Add ST_Project(geography, distance, azimuth) (#657) to construct a new point given a heading and a distance.

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@8495 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Paul Ramsey 2011-12-21 18:42:08 +00:00
parent 88a4eb1d25
commit 35c0bbf274
5 changed files with 187 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

50
doc/reference_measure.xml
View file

@ -3592,6 +3592,56 @@ SELECT ST_AsEWKT(ST_PointOnSurface(ST_GeomFromEWKT('LINESTRING(0 5 1, 0 0 1, 0 1
</refsection>
</refentry>
<refentry id="ST_Project">
<refnamediv>
<refname>ST_Project</refname>
<refpurpose>Returns a <varname>POINT</varname> projected from a start point using a bearing and distance.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<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>
</funcsynopsis>
</refsynopsisdiv>
<refsection>
<title>Description</title>
<para>Returns a <varname>POINT</varname> projected from a start point using an azimuth (bearing) and distance.</para>
<para>Distance, azimuth and projection are all aspects of the same operation, describing (or in the case of projection, constructing) the relationship between two points on the world.</para>
<para>The azimuth is sometimes called the heading or the bearing in navigation. It is measured relative to true north (azimuth zero). East is azimuth 90, south is azimuth 180, west is azimuth 270.</para>
<para>The distance is given in meters.</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting>SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, 45));
st_astext
------------------------------------------
POINT(0.63523102912532 0.63947233472882)
(1 row)
</programlisting>
</refsection>
<refsection>
<title>See Also</title>
<para><xref linkend="ST_Azimuth" />, <xref linkend="ST_Distance" /></para>
</refsection>
</refentry>
<refentry id="ST_Relate">
<refnamediv>
<refname>ST_Relate</refname>

5
liblwgeom/liblwgeom.h.in
View file

@ -1321,6 +1321,11 @@ extern void spheroid_init(SPHEROID *s, double a, double b);
*/
extern double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance);
/**
* Calculate the location of a point, give a start point, bearing and distance.
*/
extern LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth);
/**
* Calculate the geodetic area of a lwgeom on the sphere. The result
* will be multiplied by the average radius of the supplied spheroid.

63
liblwgeom/lwgeodetic.c
View file

@ -74,6 +74,7 @@ double latitude_radians_normalize(double lat)
/**
* Convert a longitude to the range of -180,180
* @param lon longitude in degrees
*/
double longitude_degrees_normalize(double lon)
{
@ -99,6 +100,7 @@ double longitude_degrees_normalize(double lon)
/**
* Convert a latitude to the range of -90,90
* @param lat latitude in degrees
*/
double latitude_degrees_normalize(double lat)
{
@ -142,6 +144,11 @@ int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *
return FP_EQUALS(g1->lat, g2->lat) && FP_EQUALS(g1->lon, g2->lon);
}
/**
* Initialize a geographic point
* @param lon longitude in degrees
* @param lat latitude in degrees
*/
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
{
g->lat = latitude_radians_normalize(deg2rad(lat));
@ -1781,6 +1788,62 @@ double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
}
/**
* Calculate a projected point given a source point, a distance and a bearing.
* @param r - location of first point.
* @param spheroid - spheroid definition.
* @param distance - distance, in units of the spheroid def'n.
* @param azimuth - azimuth in degrees.
* @return s - location of projected point.
*
*/
LWPOINT* lwgeom_project_spheroid(const LWPOINT *r, const SPHEROID *spheroid, double distance, double azimuth)
{
GEOGRAPHIC_POINT geo_source, geo_dest;
POINT4D pt_dest;
double azimuth_radians;
double x, y;
POINTARRAY *pa;
LWPOINT *lwp;
/* Check the azimuth validity, convert to radians */
if ( azimuth < -360.0 || azimuth > 360.0 )
{
lwerror("Azimuth must be between -360 and 360");
return NULL;
}
azimuth_radians = deg2rad(azimuth);
/* Check the distance validity */
if ( distance < 0.0 || distance > (M_PI * spheroid->radius) )
{
lwerror("Distance must be between 0 and %g", M_PI * spheroid->radius);
return NULL;
}
/* Convert to ta geodetic point */
x = lwpoint_get_x(r);
y = lwpoint_get_y(r);
geographic_point_init(x, y, &geo_source);
/* Try the projection */
if( spheroid_project(&geo_source, spheroid, distance, azimuth_radians, &geo_dest) == LW_FAILURE )
{
LWDEBUGF(3, "Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
lwerror("Unable to project from (%g %g) with azimuth %g and distance %g", x, y, azimuth, distance);
return NULL;
}
/* Build the output LWPOINT */
pa = ptarray_construct(0, 0, 1);
pt_dest.x = rad2deg(geo_dest.lon);
pt_dest.y = rad2deg(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);
return lwp;
}
/**
* Calculate the distance between two LWGEOMs, using the coordinates are
* longitude and latitude. Return immediately when the calulated distance drops

6
postgis/geography.sql.in.c
View file

@ -670,6 +670,12 @@ CREATE OR REPLACE FUNCTION ST_Length(text)
$$ SELECT ST_Length($1::geometry); $$
LANGUAGE 'SQL' IMMUTABLE STRICT;
-- Availability: 2.0.0
CREATE OR REPLACE FUNCTION ST_Project(geog geography, distance float8, azimuth float8)
RETURNS geography
AS 'MODULE_PATHNAME','geography_project'
LANGUAGE 'C' IMMUTABLE STRICT
COST 100;
-- Availability: 2.0.0
CREATE OR REPLACE FUNCTION ST_Perimeter(geog geography, use_spheroid boolean DEFAULT true)

63
postgis/geography_measurement.c
View file

@ -35,6 +35,7 @@ Datum geography_point_outside(PG_FUNCTION_ARGS);
Datum geography_covers(PG_FUNCTION_ARGS);
Datum geography_bestsrid(PG_FUNCTION_ARGS);
Datum geography_perimeter(PG_FUNCTION_ARGS);
Datum geography_project(PG_FUNCTION_ARGS);
/*
** geography_distance(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance, boolean use_spheroid)
@ -596,4 +597,66 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS)
}
/*
** geography_project(GSERIALIZED *g, distance, azimuth)
** returns point of projection given start point,
** azimuth (bearing) and distance
*/
PG_FUNCTION_INFO_V1(geography_project);
Datum geography_project(PG_FUNCTION_ARGS)
{
LWGEOM *lwgeom = NULL;
LWPOINT *lwp_projected;
GSERIALIZED *g = NULL;
GSERIALIZED *g_out = NULL;
double azimuth, distance;
SPHEROID s;
uint32_t type;
/* Get our geometry object loaded into memory. */
g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
/* Only return for points. */
type = gserialized_get_type(g);
if ( type != POINTTYPE )
{
elog(ERROR, "ST_Project(geography) is only valid for point inputs");
PG_RETURN_NULL();
}
lwgeom = lwgeom_from_gserialized(g);
/* EMPTY things cannot be projected from */
if ( lwgeom_is_empty(lwgeom) )
{
lwgeom_free(lwgeom);
elog(ERROR, "ST_Project(geography) cannot project from an empty start point");
PG_RETURN_NULL();
}
/* Read the other parameters */
distance = PG_GETARG_FLOAT8(1);
azimuth = PG_GETARG_FLOAT8(2);
/* Initialize spheroid */
spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
/* Calculate the length */
lwp_projected = lwgeom_project_spheroid(lwgeom_as_lwpoint(lwgeom), &s, distance, azimuth);
/* Something went wrong... */
if ( lwp_projected == NULL )
{
elog(ERROR, "lwgeom_project_spheroid returned null");
PG_RETURN_NULL();
}
/* Clean up, but not all the way to the point arrays */
lwgeom_free(lwgeom);
g_out = geography_serialize(lwpoint_as_lwgeom(lwp_projected));
lwpoint_free(lwp_projected);
PG_FREE_IF_COPY(g, 0);
PG_RETURN_POINTER(g_out);
}