ST_Length_Spheroid renamed to ST_LengthSpheroid,
ST_Distance_Spheroid renamed to ST_DistanceSpheroid,
ST_Distance_Sphere renamed to ST_DistanceSphere,
- ST_Length2D_Spheroid renamed to ST_Length2DSpheroid,
- ST_Point_Inside_Circle renamed to ST_PointInsideCircle,
ST_3DLength_Spheroid deprecated (use ST_LengthSpheroid)
- #2769, ST_Mem_Size renamed to ST_MemSize
- #2565, ST_SummaryStats(tablename, rastercolumn, ...)
- Canonical output for index key types
- ST_SwapOrdinates (Sandro Santilli / Boundless)
+ - #2918, Use GeographicLib functions for geodetics (Mike Toews)
- #3074, ST_Subdivide (Paul Ramsey / CartoDB)
- #3040, KNN GiST index based centroid (<<->>) and box (<<#>>)
n-D distance operators (Sandro Santilli / Boundless)
dnl If we couldn't find pg_config, display an error
if test "x$PG_CONFIG" = "x"; then
- AC_MSG_ERROR([could not find pg_config within the current path. You may need to try re-running configure with a --with-pgconfig parameter.])
+ AC_MSG_ERROR([could not find pg_config within the current path. You may need to try re-running configure with a --with-pg_config parameter.])
fi
else
dnl PG_CONFIG was specified; display a message to the user
</question>
<answer>
- <para>You can store point, line, polygon, multipoint, multiline,
- multipolygon, and geometrycollections. In PostGIS 2.0 and above you can also store TINS and Polyhedral Surfaces in the basic geometry type.
+ <para>You can store Point, LineString, Polygon, MultiPoint, MultiLineString,
+ MultiPolygon, and GeometryCollection geometries. In PostGIS 2.0 and above you can also store TINS and Polyhedral Surfaces in the basic geometry type.
These are specified in the Open
- GIS Well Known Text Format (with XYZ,XYM,XYZM extensions). There are three data types currently supported.
+ GIS Well Known Text Format (with Z, M, and ZM extensions). There are three data types currently supported.
The standard OGC geometry data type which uses a planar coordinate system for measurement, the
- geography data type which uses a geodetic coordinate system (not OGC, but you'll find a similar type in Microsoft SQL Server 2008+). Only WGS 84 long lat (SRID:4326) is supported
- by the geography data type. The newest family member of the PostGIS spatial type family is raster for storing and analyzing raster data. Raster has its very own FAQ. Refer to <xref linkend="RT_FAQ"/>
+ geography data type which uses a geodetic coordinate system, with calculations on either a sphere or spheroid.
+ The newest family member of the PostGIS spatial type family is raster for storing and analyzing raster data. Raster has its very own FAQ. Refer to <xref linkend="RT_FAQ"/>
and <xref linkend="RT_reference" /> for more details.</para>
</answer>
</qandaentry>
</question>
<answer>
- <para>Short Answer: geography is a new data type that supports long range distances measurements, but most computations on it are currently
+ <para>Short Answer: geography is a newer data type that supports long range distances measurements, but most computations on it are
slower than they are on geometry. If
- you use geography -- you don't need to learn much about planar coordinate systems. Geography is generally best
+ you use geography, you don't need to learn much about planar coordinate systems. Geography is generally best
if all you care about is measuring distances and lengths and you have data from all over the world.
Geometry data type is an older data type that has many more functions supporting it, enjoys greater support from third party tools,
and operations on it are generally faster -- sometimes as much as 10 fold faster for larger geometries.
Maxime Guillaud,
Maxime van Noppen,
Michael Fuhr,
+Mike Toews,
Nathan Wagner,
Nathaniel Clay,
Nikita Shulga,
</refsynopsisdiv>
<refsection>
<title>Description</title>
- <para>Returns a geography object from the well-known text or extended well-known representation. SRID 4326 is assumed. This
- is an alias for ST_GeographyFromText. Points are always expressed in long lat form.</para>
+ <para>Returns a geography object from the well-known text or extended well-known representation. SRID 4326 is assumed if unspecified.
+ This is an alias for ST_GeographyFromText. Points are always expressed in long lat form.</para>
</refsection>
<refsection>
<programlisting>
--- converting lon lat coords to geography
ALTER TABLE sometable ADD COLUMN geog geography(POINT,4326);
-UPDATE sometable SET geog = ST_GeogFromText('SRID=4326;POINT(' || lon || ' ' || lat || ')');
+UPDATE sometable SET geog = ST_GeogFromText('SRID=4326;POINT(' || lon || ' ' || lat || ')');
+
+--- specify a geography point using EPSG:4267, NAD27
+SELECT ST_AsEWKT(ST_GeogFromText('SRID=4267;POINT(-77.0092 38.889588)'));
</programlisting>
</refsection>
<refsection>
</refsynopsisdiv>
<refsection>
<title>Description</title>
- <para>Returns a geography object from the well-known text representation. SRID 4326 is assumed.</para>
+ <para>Returns a geography object from the well-known text representation. SRID 4326 is assumed if unspecified.</para>
<!-- TODO: put in example -->
</refsection>
<refsection>
<refnamediv>
<refname>ST_Area</refname>
- <refpurpose>Returns the area of the surface if it is a polygon or
- multi-polygon. For "geometry" type area is in SRID units. For "geography" area is in square meters.</refpurpose>
+ <refpurpose>Returns the area of the surface if it is a Polygon or
+ MultiPolygon. For geometry, a 2D Cartesian area is determined with units specified by the SRID. For geography, area is determined on a curved surface with units in square meters.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<refsection>
<title>Description</title>
- <para>Returns the area of the geometry if it is a polygon or
- multi-polygon. Return the area measurement of an ST_Surface or
- ST_MultiSurface value. For geometry Area is in the units of the srid. For geography area is in square meters and defaults to measuring about the spheroid of the geography (currently only WGS84).
- To measure around the faster but less accurate sphere -- ST_Area(geog,false).
+ <para>Returns the area of the geometry if it is a Polygon or
+ MultiPolygon. Return the area measurement of an ST_Surface or
+ ST_MultiSurface value. For geometry, a 2D Cartesian area is determined with units specified by the SRID. For geography, by default area is determined on a spheroid with units in square meters.
+ To measure around the faster but less accurate sphere, use ST_Area(geog,false).
</para>
<para>Enhanced: 2.0.0 - support for 2D polyhedral surfaces was introduced.</para>
+ <para>Enhanced: 2.2.0 - measurement on spheroid performed with GeographicLib for improved accuracy and robustness.</para>
<para>&sfs_compliant;</para>
<para>&sqlmm_compliant; SQL-MM 3: 8.1.2, 9.5.3</para>
<para>&P_support;</para>
<refsection>
<title>Examples</title>
<para>Return area in square feet for a plot of Massachusetts land and multiply by conversion to get square meters.
- Note this is in square feet because 2249 is
- Mass State Plane Feet </para>
+ Note this is in square feet because EPSG:2249 is
+ Massachusetts State Plane Feet </para>
<programlisting>
SELECT ST_Area(the_geom) As sqft, ST_Area(the_geom)*POWER(0.3048,2) As sqm
FROM (SELECT
---------+-------------
928.625 | 86.27208552
</programlisting>
-<para>Return area square feet and transform to Massachusetts state plane meters (26986) to get square meters.
+<para>Return area square feet and transform to Massachusetts state plane meters (EPSG:26986) to get square meters.
Note this is in square feet because 2249 is
- Mass State Plane Feet and transformed area is in square meters since 26986 is state plane mass meters </para>
+ Massachusetts State Plane Feet and transformed area is in square meters since EPSG:26986 is state plane Massachusetts meters </para>
<programlisting>
SELECT ST_Area(the_geom) As sqft, ST_Area(ST_Transform(the_geom,26986)) As sqm
928.625 | 86.2724304199219
</programlisting>
-<para>Return area square feet and square meters using Geography data type. Note that we transform to our geometry to geography
+<para>Return area square feet and square meters using geography data type. Note that we transform to our geometry to geography
(before you can do that make sure your geometry is in WGS 84 long lat 4326). Geography always measures in meters.
This is just for demonstration to compare. Normally your table will be stored in geography data type already.</para>
<programlisting>
)
)
) As foo(the_geog);
- sqft_spheroid | sqft_sphere | sqm_spheroid
------------------+------------------+------------------
-928.684405217197 | 927.186481558724 | 86.2776044452694
+ sqft_spheroid | sqft_sphere | sqm_spheroid
+------------------+------------------+------------------
+ 928.684403538925 | 927.049336105925 | 86.2776042893529
--if your data is in geography already
SELECT ST_Area(the_geog)/POWER(0.3048,2) As sqft, ST_Area(the_geog) As sqm
<title>Description</title>
<para>Returns the azimuth in radians of the segment defined by the given
- point-geometries, or NULL if the two points are coincident. The azimuth is north-based and is measured clockwise: North = 0; East = PI/2; South = PI; West = 3PI/2.</para>
-
- <para>The Azimuth is mathematical concept defined as the angle, in this case measured in radian, between a reference plane
- and a point. </para>
+ point geometries, or NULL if the two points are coincident. The azimuth is angle is referenced from north, and is positive clockwise: North = 0; East = π/2; South = π; West = 3π/2.</para>
+ <para>For the geography type, the forward azimuth is solved as part of the inverse geodesic problem.</para>
+ <para>The azimuth is mathematical concept defined as the angle between a reference plane and a point, with angular units in radians.
+ Units can be converted to degrees using a built-in PostgreSQL function degrees(), as shown in the example.</para>
<para>Availability: 1.1.0</para>
<para>Enhanced: 2.0.0 support for geography was introduced.</para>
+ <para>Enhanced: 2.2.0 measurement on spheroid performed with GeographicLib for improved accuracy and robustness.</para>
<para>Azimuth is especially useful in conjunction with ST_Translate for shifting an object along its perpendicular axis. See
upgis_lineshift <ulink url="http://trac.osgeo.org/postgis/wiki/UsersWikiplpgsqlfunctions">Plpgsqlfunctions PostGIS wiki section</ulink> for example of this.</para>
</refsection>
<title>Examples</title>
<para>Geometry Azimuth in degrees </para>
<programlisting>
-SELECT ST_Azimuth(ST_Point(25,45), ST_Point(75,100))/(2*pi())*360 as degA_B,
- ST_Azimuth(ST_Point(75,100), ST_Point(25,45))/(2*pi())*360 As degB_A;
-
--- NOTE easier to remember syntax using PostgreSQL built-in degrees function --
--- Both yield same answer --
-SELECT degrees( ST_Azimuth(ST_Point(25,45), ST_Point(75,100)) ) as degA_B,
- degrees( ST_Azimuth(ST_Point(75,100), ST_Point(25,45)) ) As degB_A;
+SELECT degrees(ST_Azimuth(ST_Point(25, 45), ST_Point(75, 100))) AS degA_B,
+ degrees(ST_Azimuth(ST_Point(75, 100), ST_Point(25, 45))) AS degB_A;
dega_b | degb_a
------------------+------------------
<refnamediv>
<refname>ST_Distance</refname>
- <refpurpose>For geometry type Returns the 2-dimensional cartesian minimum distance (based on spatial ref) between two geometries in
- projected units. For geography type defaults to return spheroidal minimum distance between two geographies in meters.</refpurpose>
+ <refpurpose>For geometry type Returns the 2D Cartesian distance between two geometries in
+ projected units (based on spatial ref). For geography type defaults to return minimum geodesic distance between two geographies in meters.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<refsection>
<title>Description</title>
- <para>For geometry type returns the 2-dimensional minimum cartesian distance between two geometries in
- projected units (spatial ref units). For geography type defaults to return the minimum distance around WGS 84 spheroid between two geographies in meters. Pass in
- false to return answer in sphere instead of spheroid.</para>
+ <para>For geometry type returns the minimum 2D Cartesian distance between two geometries in
+ projected units (spatial ref units). For geography type defaults to return the minimum geodesic distance between two geographies in meters. If use_spheroid is
+ false, a faster sphere calculation is used instead of a spheroid.</para>
<para>&sfs_compliant;</para>
<para>&sqlmm_compliant; SQL-MM 3: 5.1.23</para>
<para>&curve_support;</para>
-
<para>&sfcgal_enhanced;</para>
+ <para>&sfcgal_enhanced;</para>
<para>Availability: 1.5.0 geography support was introduced in 1.5. Speed improvements for planar to better handle large or many vertex geometries</para>
- <para>Enhanced: 2.1.0 improved speed for geography. See <ulink url="http://blog.opengeo.org/2012/07/12/making-geography-faster/">Making Geography faster</ulink> for details.</para>
+ <para>Enhanced: 2.1.0 improved speed for geography. See <ulink url="http://boundlessgeo.com/2012/07/making-geography-faster/">Making Geography faster</ulink> for details.</para>
<para>Enhanced: 2.1.0 - support for curved geometries was introduced.</para>
+ <para>Enhanced: 2.2.0 - measurement on spheroid performed with GeographicLib for improved accuracy and robustness.</para>
</refsection>
<refsection>
having the same SRID.</para>
<para>For geography units are in meters and measurement is
- defaulted to use_spheroid=true (measure around WGS 84 spheroid), for faster check, use_spheroid=false to measure along sphere.
+ defaulted to use_spheroid=true, for faster check, use_spheroid=false to measure along sphere.
</para>
<note>
<para>This function call will automatically include a bounding box
<refnamediv>
<refname>ST_Length</refname>
- <refpurpose>Returns the 2d length of the geometry if it is a linestring or multilinestring. geometry are in units of spatial reference and geography are in meters (default spheroid)</refpurpose>
+ <refpurpose>Returns the 2D length of the geometry if it is a LineString or MultiLineString. geometry are in units of spatial reference and geography are in meters (default spheroid)</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<refsection>
<title>Description</title>
- <para>For geometry: Returns the cartesian 2D length of the geometry if it is a linestring, multilinestring, ST_Curve, ST_MultiCurve. 0 is returned for
- areal geometries. For areal geometries use ST_Perimeter. Geometry: Measurements are in the units of the
- spatial reference system of the geometry. Geography: Units are in meters and also acts as a Perimeter function for areal geogs.</para>
+ <para>For geometry: Returns the 2D Cartesian length of the geometry if it is a LineString, MultiLineString, ST_Curve, ST_MultiCurve. 0 is returned for
+ areal geometries. For areal geometries use <xref linkend="ST_Perimeter" />. For geometry types, units for length measures are specified by the
+ spatial reference system of the geometry.</para>
+ <para>For geography types, the calculations are performed using the inverse geodesic problem, where length units are in meters.
+ If PostGIS is compiled with PROJ version 4.8.0 or later, the spheroid is specified by the SRID, otherwise it is exclusive to WGS84.
+ If <varname>use_spheroid=false</varname>, then calculations will approximate a sphere instead of a spheroid.</para>
<para>Currently for geometry this is an alias for ST_Length2D, but this may change to support higher dimensions.</para>
<warning><para>Changed: 2.0.0 Breaking change -- in prior versions applying this to a MULTI/POLYGON of type geography would give you the perimeter of the POLYGON/MULTIPOLYGON. In 2.0.0
<para>&sfs_compliant; s2.1.5.1</para>
<para>&sqlmm_compliant; SQL-MM 3: 7.1.2, 9.3.4</para>
<para>Availability: 1.5.0 geography support was introduced in 1.5.</para>
-
<para>&sfcgal_enhanced;</para>
+ <para>&sfcgal_enhanced;</para>
</refsection>
<refsection>
<title>Geometry Examples</title>
- <para>Return length in feet for line string. Note this is in feet because 2249 is
- Mass State Plane Feet</para>
+ <para>Return length in feet for line string. Note this is in feet because EPSG:2249 is
+ Massachusetts State Plane Feet</para>
<programlisting>
SELECT ST_Length(ST_GeomFromText('LINESTRING(743238 2967416,743238 2967450,743265 2967450,
743265.625 2967416,743238 2967416)',2249));
122.630744000095
---Transforming WGS 84 linestring to Massachusetts state plane meters
+--Transforming WGS 84 LineString to Massachusetts state plane meters
SELECT ST_Length(
ST_Transform(
ST_GeomFromEWKT('SRID=4326;LINESTRING(-72.1260 42.45, -72.1240 42.45666, -72.123 42.1546)'),
As foo;
length_spheroid | length_sphere
------------------+------------------
- 34310.5703627305 | 34346.2060960742
-(1 row)
+ 34310.5703627288 | 34346.2060960742
</programlisting>
</refsection>
<refsection>
<refsection>
<title>Examples</title>
- <para>Return length in feet for a 3D cable. Note this is in feet because 2249 is
- Mass State Plane Feet</para>
+ <para>Return length in feet for a 3D cable. Note this is in feet because EPSG:2249 is
+ Massachusetts State Plane Feet</para>
<programlisting>
SELECT ST_3DLength(ST_GeomFromText('LINESTRING(743238 2967416 1,743238 2967450 1,743265 2967450 3,
743265.625 2967416 3,743238 2967416 3)',2249));
<refentry id="ST_Length2D_Spheroid">
<refnamediv>
- <refname>ST_Length2DSpheroid</refname>
+ <refname>ST_Length2D_Spheroid</refname>
<refpurpose>Calculates the 2D length of a linestring/multilinestring on an ellipsoid. This
is useful if the coordinates of the geometry are in
<refsynopsisdiv>
<funcsynopsis>
<funcprototype>
- <funcdef>float <function>ST_Length2DSpheroid</function></funcdef>
+ <funcdef>float <function>ST_Length2D_Spheroid</function></funcdef>
<paramdef><type>geometry </type> <parameter>a_linestring</parameter></paramdef>
<paramdef><type>spheroid </type> <parameter>a_spheroid</parameter></paramdef>
</funcprototype>
<note><para>Will return 0 for anything that is not a MULTILINESTRING or LINESTRING</para></note>
<note><para>This is much like <xref linkend="ST_Length_Spheroid" /> except it will throw away the Z coordinate in calculations.</para></note>
- <para>Availability: 1.2.2</para>
- <para>Changed: 2.2.0 In prior versions this used to be called ST_Length2D_Spheroid</para>
-
</refsection>
<refsection>
<title>Examples</title>
- <programlisting>SELECT ST_Length2DSpheroid( geometry_column,
+ <programlisting>SELECT ST_Length2D_Spheroid( geometry_column,
'SPHEROID["GRS_1980",6378137,298.257222101]' )
FROM geometry_table;
-SELECT ST_Length2DSpheroid( the_geom, sph_m ) As tot_len,
-ST_Length2DSpheroid(ST_GeometryN(the_geom,1), sph_m) As len_line1,
-ST_Length2DSpheroid(ST_GeometryN(the_geom,2), sph_m) As len_line2
+SELECT ST_Length2D_Spheroid( the_geom, sph_m ) As tot_len,
+ST_Length2D_Spheroid(ST_GeometryN(the_geom,1), sph_m) As len_line1,
+ST_Length2D_Spheroid(ST_GeometryN(the_geom,2), sph_m) As len_line2
FROM (SELECT ST_GeomFromText('MULTILINESTRING((-118.584 38.374,-118.583 38.5),
(-71.05957 42.3589 , -71.061 43))') As the_geom,
CAST('SPHEROID["GRS_1980",6378137,298.257222101]' As spheroid) As sph_m) as foo;
85204.5207562955 | 13986.8725229309 | 71217.6482333646
--3D Observe same answer
-SELECT ST_Length2DSpheroid( the_geom, sph_m ) As tot_len,
-ST_Length2DSpheroid(ST_GeometryN(the_geom,1), sph_m) As len_line1,
-ST_Length2DSpheroid(ST_GeometryN(the_geom,2), sph_m) As len_line2
+SELECT ST_Length2D_Spheroid( the_geom, sph_m ) As tot_len,
+ST_Length2D_Spheroid(ST_GeometryN(the_geom,1), sph_m) As len_line1,
+ST_Length2D_Spheroid(ST_GeometryN(the_geom,2), sph_m) As len_line2
FROM (SELECT ST_GeomFromEWKT('MULTILINESTRING((-118.584 38.374 20,-118.583 38.5 30),
(-71.05957 42.3589 75, -71.061 43 90))') As the_geom,
CAST('SPHEROID["GRS_1980",6378137,298.257222101]' As spheroid) As sph_m) as foo;
<refname>ST_Perimeter</refname>
<refpurpose>Return the length measurement of the boundary of an ST_Surface
- or ST_MultiSurface geometry or geography. (Polygon, Multipolygon). geometry measurement is in units of spatial reference and geography is in meters.</refpurpose>
+ or ST_MultiSurface geometry or geography. (Polygon, MultiPolygon). geometry measurement is in units of spatial reference and geography is in meters.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<refsection>
<title>Description</title>
- <para>Returns the 2D perimeter of the geometry/geography if it is a ST_Surface, ST_MultiSurface (Polygon, Multipolygon). 0 is returned for
- non-areal geometries. For linestrings use ST_Length. Measurements for geometry are in the units of the
- spatial reference system of the geometry. Measurements for geography are in meters. If <varname>use_spheroid</varname> is set to false, then will
- model earth as a sphere instead of a spheroid.</para>
+ <para>Returns the 2D perimeter of the geometry/geography if it is a ST_Surface, ST_MultiSurface (Polygon, MultiPolygon). 0 is returned for
+ non-areal geometries. For linear geometries use <xref linkend="ST_Length" />. For geometry types, units for perimeter measures are specified by the
+ spatial reference system of the geometry.</para>
+ <para>For geography types, the calculations are performed using the inverse geodesic problem, where perimeter units are in meters.
+ If PostGIS is compiled with PROJ version 4.8.0 or later, the spheroid is specified by the SRID, otherwise it is exclusive to WGS84.
+ If <varname>use_spheroid=false</varname>, then calculations will approximate a sphere instead of a spheroid.</para>
<para>Currently this is an alias for ST_Perimeter2D, but this may change to support higher dimensions.</para>
<refsection>
<title>Examples: Geometry</title>
- <para>Return perimeter in feet for polygon and multipolygon. Note this is in feet because 2249 is
- Mass State Plane Feet</para>
+ <para>Return perimeter in feet for Polygon and MultiPolygon. Note this is in feet because EPSG:2249 is
+ Massachusetts State Plane Feet</para>
<programlisting>
SELECT ST_Perimeter(ST_GeomFromText('POLYGON((743238 2967416,743238 2967450,743265 2967450,
743265.625 2967416,743238 2967416))', 2249));
</refsection>
<refsection>
<title>Examples: Geography</title>
- <para>Return perimeter in meters and feet for polygon and multipolygon. Note this is geography (WGS 84 long lat)</para>
+ <para>Return perimeter in meters and feet for Polygon and MultiPolygon. Note this is geography (WGS 84 long lat)</para>
<programlisting>
SELECT ST_Perimeter(geog) As per_meters, ST_Perimeter(geog)/0.3048 As per_ft
FROM ST_GeogFromText('POLYGON((-71.1776848522251 42.3902896512902,-71.1776843766326 42.3903829478009,
37.3790462565251 | 122.634666195949
--- Multipolygon example --
+-- MultiPolygon example --
SELECT ST_Perimeter(geog) As per_meters, ST_Perimeter(geog,false) As per_sphere_meters, ST_Perimeter(geog)/0.3048 As per_ft
FROM ST_GeogFromText('MULTIPOLYGON(((-71.1044543107478 42.340674480411,-71.1044542869917 42.3406744369506,
-71.1044553562977 42.340673886454,-71.1044543107478 42.340674480411)),
<refsection>
<title>Description</title>
- <para>Returns a <varname>POINT</varname> projected from a start point using an azimuth (bearing) measured in radians and distance measured in meters.</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 (<varname>pi/2</varname>), south is azimuth 180 (<varname>pi</varname>), west is azimuth 270 (<varname>pi*1.5</varname>).</para>
+ <para>Returns a <varname>POINT</varname> projected from a start point using an azimuth (bearing) measured in radians and distance measured in meters. This is also called a direct geodesic problem.</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 (π/2), south is azimuth 180 (π), west is azimuth 270 (3π/2).</para>
<para>The distance is given in meters.</para>
<para>Availability: 2.0.0</para>
<programlisting>SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, radians(45.0)));
st_astext
-------------------------------------------
- POINT(0.63523102912532 0.63947233472882)
+ st_astext
+--------------------------------------------
+ POINT(0.635231029125537 0.639472334729198)
(1 row)
</programlisting>
</refsection>
<refsection>
- <title>Example: Using radians - projected point 100,000 meters and bearing pi/4 (45 degrees) </title>
+ <title>Example: Using radians - projected point 100,000 meters and bearing pi/4 radians (45 degrees) </title>
<programlisting>SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, pi()/4));
- st_astext
-------------------------------------------
- POINT(0.63523102912532 0.63947233472882)
+ st_astext
+--------------------------------------------
+ POINT(0.635231029125537 0.639472334729198)
(1 row)
</programlisting>
</refsection>
If g1 and g2 are intersecting with more than one point the function will return a line with start
and end in the same point but it can be any of the intersecting points.
The line returned will always start in g1 and end in g2.
- The length of the line this function returns will always be the same as st_distance returns for g1 and g2.
+ The length of the line this function returns will always be the same as ST_Distance returns for g1 and g2.
</para>
<para>Availability: 1.5.0</para>
varint.o
NM_OBJS = \
+ geodesic.o \
lwspheroid.o
ifeq (@SFCGAL@,sfcgal)
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
*
+ * Note: Geodesic measurements have been independently verified using
+ * GeographicLib v 1.37 with MPFR C++ using utilities GeodSolve and
+ * Planimeter, with -E "exact" flag with the following env vars:
+ * export GEOGRAPHICLIB_DIGITS=1000
+ * export WGS84_ELIPSOID="6378137 298.257223563"
+ * export WGS84_SPHERE="6371008.7714150598325213222 0"
**********************************************************************/
#include <stdio.h>
geographic_point_init(1, 0, &e);
dist = sphere_distance(&s, &e);
dir = sphere_direction(&s, &e, dist);
- CU_ASSERT_DOUBLE_EQUAL(dir, M_PI / 2.0, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "0 0 0 1" */
+ CU_ASSERT_DOUBLE_EQUAL(dir, M_PI_2, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist, 0.0174532925199433, 1e-14);
geographic_point_init(0, 0, &s);
geographic_point_init(0, 1, &e);
dist = sphere_distance(&s, &e);
dir = sphere_direction(&s, &e, dist);
- CU_ASSERT_DOUBLE_EQUAL(dir, 0.0, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "0 0 1 0" */
+ CU_ASSERT_DOUBLE_EQUAL(dir, 0.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist, 0.0174532925199433, 1e-14);
}
GEOGRAPHIC_POINT s, e;
double dir1, dist1, dir2, dist2;
- dir1 = M_PI/2;
+ dir1 = M_PI_2;
dist1 = 0.1;
- geographic_point_init(0, 0, &s);
+ geographic_point_init(0, 0, &s);
sphere_project(&s, dist1, dir1, &e);
+
+ CU_ASSERT_DOUBLE_EQUAL(e.lon, 0.1, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(e.lat, 0.0, 1e-14);
+ /* Direct and inverse solutions agree */
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist1);
- CU_ASSERT_DOUBLE_EQUAL(dist1, dist2, 0.0001);
- CU_ASSERT_DOUBLE_EQUAL(dir1, dir2, 0.0001);
+ CU_ASSERT_DOUBLE_EQUAL(dist1, dist2, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dir1, dir2, 1e-14);
dist1 = sphere_distance(&e, &s);
dir1 = sphere_direction(&e, &s, dist1);
sphere_project(&e, dist1, dir1, &s);
-
- CU_ASSERT_DOUBLE_EQUAL(s.lon, 0.0, 0.0001);
- CU_ASSERT_DOUBLE_EQUAL(s.lat, 0.0, 0.0001);
+
+ CU_ASSERT_DOUBLE_EQUAL(s.lon, 0.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(s.lat, 0.0, 1e-14);
geographic_point_init(0, 0.2, &e);
geographic_point_init(0, 0.4, &s);
dist1 = sphere_distance(&s, &e);
dir1 = sphere_direction(&e, &s, dist1);
- CU_ASSERT_DOUBLE_EQUAL(dir1, 0.0, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "0.2 0 0.4 0" */
+ CU_ASSERT_DOUBLE_EQUAL(dir1, 0.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist1, 0.0034906585039887, 1e-14);
- geographic_point_init(0, 1, &s);
+ geographic_point_init(0, 1, &s); /* same start point for remainder of tests */
geographic_point_init(0, 2, &e);
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist2);
- CU_ASSERT_DOUBLE_EQUAL(dir2, 0.0, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 2 0" */
+ CU_ASSERT_DOUBLE_EQUAL(dir2, 0.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174532925199433, 1e-14);
geographic_point_init(1, 1, &e);
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist2);
- CU_ASSERT_DOUBLE_EQUAL(dir2, 1.57064, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 1 1" */
+ CU_ASSERT_DOUBLE_EQUAL(dir2, 89.991273575329292895136 * M_PI / 180.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174506342314906, 1e-14);
geographic_point_init(0, 0, &e);
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist2);
- CU_ASSERT_DOUBLE_EQUAL(dir2, 3.14159, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 0 0" */
+ CU_ASSERT_DOUBLE_EQUAL(dir2, M_PI, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174532925199433, 1e-14);
geographic_point_init(-1, 1, &e);
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist2);
- CU_ASSERT_DOUBLE_EQUAL(dir2, -1.57064, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 1 -1" */
+ CU_ASSERT_DOUBLE_EQUAL(dir2, -89.991273575329292895136 * M_PI / 180.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0174506342314906, 1e-14);
geographic_point_init(1, 2, &e);
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist2);
- CU_ASSERT_DOUBLE_EQUAL(dir2, 0.785017, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 2 1" */
+ CU_ASSERT_DOUBLE_EQUAL(dir2, 44.978182941465044354783 * M_PI / 180.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0246782972905467, 1e-14);
geographic_point_init(-1, 0, &e);
dist2 = sphere_distance(&s, &e);
dir2 = sphere_direction(&s, &e, dist2);
- CU_ASSERT_DOUBLE_EQUAL(dir2, -2.35612, 0.0001);
+ /* GeodSolve -i -E -p 16 -e 1 0 --input-string "1 0 0 -1" */
+ CU_ASSERT_DOUBLE_EQUAL(dir2, -134.995636455344851488216 * M_PI / 180.0, 1e-14);
+ CU_ASSERT_DOUBLE_EQUAL(dist2, 0.0246820563917664, 1e-14);
}
#if 0
gbox_geocentric_slow = LW_FALSE;
if (
- ( fabs( gbox.xmin - gbox_slow.xmin ) > gtolerance ) ||
- ( fabs( gbox.xmax - gbox_slow.xmax ) > gtolerance ) ||
- ( fabs( gbox.ymin - gbox_slow.ymin ) > gtolerance ) ||
- ( fabs( gbox.ymax - gbox_slow.ymax ) > gtolerance ) ||
- ( fabs( gbox.zmin - gbox_slow.zmin ) > gtolerance ) ||
- ( fabs( gbox.zmax - gbox_slow.zmax ) > gtolerance ) )
+ ( fabs( gbox.xmin - gbox_slow.xmin ) > gtolerance ) ||
+ ( fabs( gbox.xmax - gbox_slow.xmax ) > gtolerance ) ||
+ ( fabs( gbox.ymin - gbox_slow.ymin ) > gtolerance ) ||
+ ( fabs( gbox.ymax - gbox_slow.ymax ) > gtolerance ) ||
+ ( fabs( gbox.zmin - gbox_slow.zmin ) > gtolerance ) ||
+ ( fabs( gbox.zmax - gbox_slow.zmax ) > gtolerance ) )
{
printf("\n-------\n");
printf("If you are seeing this, cut and paste it, it is a randomly generated test case!\n");
CU_ASSERT_DOUBLE_EQUAL(closest.lon, 0.0, 0.00001);
/* Ticket #2351 */
- edge_set(149.386990599235, -26.3567415843982, 149.386990599247, -26.3567415843965, &e);
+ edge_set(149.386990599235, -26.3567415843982, 149.386990599247, -26.3567415843965, &e);
point_set(149.386990599235, -26.3567415843982, &g);
d = edge_distance_to_point(&e, &g, &closest);
CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
- // printf("CLOSE POINT(%g %g)\n", closest.lon, closest.lat);
- // printf(" ORIG POINT(%g %g)\n", g.lon, g.lat);
+ // printf("CLOSE POINT(%g %g)\n", closest.lon, closest.lat);
+ // printf(" ORIG POINT(%g %g)\n", g.lon, g.lat);
CU_ASSERT_DOUBLE_EQUAL(g.lat, closest.lat, 0.00001);
CU_ASSERT_DOUBLE_EQUAL(g.lon, closest.lon, 0.00001);
}
result = lwpoly_covers_point2d(poly, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
+
+ /* Triangle over the antimeridian */
+ lwg = lwgeom_from_wkt("POLYGON((140 52, 152.0 -6.0, -120.0 -29.0, 140 52))", LW_PARSER_CHECK_NONE);
+ poly = (LWPOLY*)lwg;
+ pt_to_test.x = -172.0;
+ pt_to_test.y = -13.0;
+ result = lwpoly_covers_point2d(poly, &pt_to_test);
+ CU_ASSERT_EQUAL(result, LW_TRUE);
+ lwgeom_free(lwg);
}
lwgeom_free(lwg1);
lwgeom_free(lwg2);
- /* Ticket #2351 */
+ /* Ticket #2351 */
lwg1 = lwgeom_from_wkt("LINESTRING(149.386990599235 -26.3567415843982,149.386990599247 -26.3567415843965)", LW_PARSER_CHECK_NONE);
lwg2 = lwgeom_from_wkt("POINT(149.386990599235 -26.3567415843982)", LW_PARSER_CHECK_NONE);
d = lwgeom_distance_spheroid(lwg1, lwg2, &s, 0.0);
{
GEOGRAPHIC_POINT g1, g2;
double d;
+#ifdef USE_PRE22GEODESIC
+ double epsilon; /* irregular */
+#else
+ const double epsilon = 1e-8; /* at least 10 nm precision */
+#endif
SPHEROID s;
/* Init to WGS84 */
spheroid_init(&s, 6378137.0, 6356752.314245179498);
- /* One vertical degree */
+ /* One vertical degree
+ $ GeodSolve -E -i -p 16 --input-string "0 0 1 0" */
point_set(0.0, 0.0, &g1);
point_set(0.0, 1.0, &g2);
d = spheroid_distance(&g1, &g2, &s);
- CU_ASSERT_DOUBLE_EQUAL(d, 110574.388615329, 0.001);
+#ifdef USE_PRE22GEODESIC
+ epsilon = 1e-6;
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(d, 110574.3885577987957342, epsilon);
- /* Ten horizontal degrees */
+ /* Ten horizontal degrees
+ $ GeodSolve -E -i -p 16 --input-string "0 -10 0 0" */
point_set(-10.0, 0.0, &g1);
point_set(0.0, 0.0, &g2);
d = spheroid_distance(&g1, &g2, &s);
- CU_ASSERT_DOUBLE_EQUAL(d, 1113194.90793274, 0.001);
+#ifdef USE_PRE22GEODESIC
+ epsilon = 1e-3;
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(d, 1113194.9079327357264771, epsilon);
- /* One horizonal degree */
+ /* One horizonal degree
+ $ GeodSolve -E -i -p 16 --input-string "0 -1 0 0" */
point_set(-1.0, 0.0, &g1);
point_set(0.0, 0.0, &g2);
d = spheroid_distance(&g1, &g2, &s);
- CU_ASSERT_DOUBLE_EQUAL(d, 111319.490779, 0.001);
+#ifdef USE_PRE22GEODESIC
+ epsilon = 1e-4;
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(d, 111319.4907932735726477, epsilon);
- /* Around world w/ slight bend */
+ /* Around world w/ slight bend
+ $ GeodSolve -E -i -p 16 --input-string "0 -180 1 0" */
point_set(-180.0, 0.0, &g1);
point_set(0.0, 1.0, &g2);
d = spheroid_distance(&g1, &g2, &s);
- CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0704483, 0.001);
+#ifdef USE_PRE22GEODESIC
+ epsilon = 1e-5;
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0700676468277450, epsilon);
- /* Up to pole */
+ /* Up to pole
+ $ GeodSolve -E -i -p 16 --input-string "0 -180 90 0" */
point_set(-180.0, 0.0, &g1);
point_set(0.0, 90.0, &g2);
d = spheroid_distance(&g1, &g2, &s);
- CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001);
+#ifdef USE_PRE22GEODESIC
+ epsilon = 1e-6;
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7293127228117396, epsilon);
}
/* Medford lot test polygon */
lwg = lwgeom_from_wkt("POLYGON((-122.848227067007 42.5007249610493,-122.848309475585 42.5007179884263,-122.848327688675 42.500835880696,-122.848245279942 42.5008428533324,-122.848227067007 42.5007249610493))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
+ /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string \
+ "42.5007249610493 -122.848227067007;42.5007179884263 -122.848309475585;"\
+ "42.500835880696 -122.848327688675;42.5008428533324 -122.848245279942" */
a1 = lwgeom_area_sphere(lwg, &s);
+ CU_ASSERT_DOUBLE_EQUAL(a1, 89.721147136698008, 0.1);
+ /* spheroid: Planimeter -E -p 20 -r --input-string \
+ "42.5007249610493 -122.848227067007;42.5007179884263 -122.848309475585;"\
+ "42.500835880696 -122.848327688675;42.5008428533324 -122.848245279942" */
a2 = lwgeom_area_spheroid(lwg, &s);
- //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
- CU_ASSERT_DOUBLE_EQUAL(a1, 89.7127703297, 0.1); /* sphere */
- CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.1); /* spheroid */
+ CU_ASSERT_DOUBLE_EQUAL(a2, 89.868413479309585, 0.1);
lwgeom_free(lwg);
/* Big-ass polygon */
lwg = lwgeom_from_wkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
+ /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string "3 -2;4 -2;4 -1;3 -1" */
a1 = lwgeom_area_sphere(lwg, &s);
+ CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.106982993974659, 0.1);
+ /* spheroid: Planimeter -E -p 20 -r --input-string "3 -2;4 -2;4 -1;3 -1" */
a2 = lwgeom_area_spheroid(lwg, &s);
- //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
- CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */
- CU_ASSERT_DOUBLE_EQUAL(a2, 12286574431.9, 10.0); /* spheroid */
+ CU_ASSERT_DOUBLE_EQUAL(a2, 12286884908.946891319597874, 0.1);
lwgeom_free(lwg);
/* One-degree square */
lwg = lwgeom_from_wkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
+ /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE --input-string "2 8.5;1 8.5;1 9.5;2 9.5" */
a1 = lwgeom_area_sphere(lwg, &s);
+ CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1);
+ /* spheroid: Planimeter -E -p 20 --input-string "2 8.5;1 8.5;1 9.5;2 9.5" */
a2 = lwgeom_area_spheroid(lwg, &s);
- //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
- CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
- CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
+ CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1);
lwgeom_free(lwg);
- /* One-degree square *near* dateline */
+ /* One-degree square *near* the antimeridian */
lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,178.5 1,178.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
+ /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE -r --input-string "2 179.5;1 179.5;1 178.5;2 178.5" */
a1 = lwgeom_area_sphere(lwg, &s);
+ CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1);
+ /* spheroid: Planimeter -E -p 20 -r --input-string "2 179.5;1 179.5;1 178.5;2 178.5" */
a2 = lwgeom_area_spheroid(lwg, &s);
- //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
- CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
- CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
+ CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1);
lwgeom_free(lwg);
- /* One-degree square *across* dateline */
+ /* One-degree square *across* the antimeridian */
lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
+ /* sphere: Planimeter -E -p 20 -e $WGS84_SPHERE --input-string "2 179.5;1 179.5;1 -179.5;2 -179.5" */
a1 = lwgeom_area_sphere(lwg, &s);
+ CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.368023059138681, 0.1);
+ /* spheroid: Planimeter -E -p 20 --input-string "2 179.5;1 179.5;1 -179.5;2 -179.5" */
a2 = lwgeom_area_spheroid(lwg, &s);
- //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
- CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.3679, 10.0); /* sphere */
- CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
+ CU_ASSERT_DOUBLE_EQUAL(a2, 12305128751.042900673161556, 0.1);
lwgeom_free(lwg);
}
CU_ASSERT_DOUBLE_EQUAL(a2, 0.017764, 0.0000001);
lwgeom_free(lwg);
- /* One-degree square *across* dateline */
+ /* One-degree square *across* antimeridian */
lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
a1 = gbox_angular_width(&gbox);
CU_ASSERT_DOUBLE_EQUAL(a2, 0.0174553, 0.0000001);
lwgeom_free(lwg);
- /* One-degree square *across* dateline */
+ /* One-degree square *across* antimeridian */
lwg = lwgeom_from_wkt("POLYGON((178.5 2,178.5 1,-179.5 1,-179.5 2,178.5 2))", LW_PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
gbox_centroid(&gbox, &pt);
p1.x = 1.0;
p2.y = 1.0;
angle = vector_angle(&p1, &p2);
- CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/2, 0.00001);
+ CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_2, 0.00001);
p1.x = p2.y = 0.0;
p1.y = 1.0;
p2.x = 1.0;
angle = vector_angle(&p1, &p2);
- CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/2, 0.00001);
+ CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_2, 0.00001);
p2.y = p2.x = 1.0;
normalize(&p2);
angle = vector_angle(&p1, &p2);
- CU_ASSERT_DOUBLE_EQUAL(angle, M_PI/4, 0.00001);
+ CU_ASSERT_DOUBLE_EQUAL(angle, M_PI_4, 0.00001);
p2.x = p2.y = p2.z = 1.0;
normalize(&p2);
p1.x = 1.0;
p2.y = 1.0;
- angle = M_PI/4;
+ angle = M_PI_4;
vector_rotate(&p1, &p2, angle, &n);
//printf("%g %g %g\n\n", n.x, n.y, n.z);
CU_ASSERT_DOUBLE_EQUAL(n.x, 0.707107, 0.00001);
/* Arc to right of the unit semicircle */
d = lw_arc_length(&A1, &A2, &A3);
- CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001);
+ CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
d = lw_arc_length(&A3, &A2, &A1);
- CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI/2, 0.000001);
+ CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
}
static void
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
/* Point 45 degrees off arc, 2 radii from center */
- P.y = P.x = 2 * cos(M_PI/4);
+ P.y = P.x = 2 * cos(M_PI_4);
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
CU_ASSERT_EQUAL( rv, LW_SUCCESS );
/* Initialize */
gbox->xmin = gbox->ymin = gbox->zmin = gbox->mmin = FLT_MAX;
- gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = -1 * FLT_MAX;
+ gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = FLT_MIN;
for ( i = 2; i < curve->points->npoints; i += 2 )
{
#include "liblwgeom.h"
+/* Define to enable pre-version 2.2 geodesic functions for geography types
+ (e.g. Vincenty for ST_Distance); otherwise use GeographicLib */
+#define USE_PRE22GEODESIC
+#undef USE_PRE22GEODESIC
/**
-* PI
-*/
-#define PI 3.1415926535897932384626433832795
-
-
-/**
-* Floating point comparitors.
+* Floating point comparators.
*/
#define FP_TOLERANCE 1e-12
#define FP_IS_ZERO(A) (fabs(A) <= FP_TOLERANCE)
/**
* Conversion functions
*/
-#define deg2rad(d) (PI * (d) / 180.0)
-#define rad2deg(r) (180.0 * (r) / PI)
+#define deg2rad(d) (M_PI * (d) / 180.0)
+#define rad2deg(r) (180.0 * (r) / M_PI)
/**
* Ape a java function
* only primitive objects should get a tree
*/
+
LWCOLLECTION *
lwgeom_subdivide(const LWGEOM *geom, int maxvertices)
{
+ int n = 0;
LWCOLLECTION *col;
col = lwcollection_construct_empty(COLLECTIONTYPE, geom->srid, lwgeom_has_z(geom), lwgeom_has_m(geom));
- lwgeom_subdivide_recursive(geom, maxvertices, col, NULL);
+ n = lwgeom_subdivide_recursive(geom, maxvertices, col, NULL);
return col;
}
#include "lwgeodetic.h"
#include "lwgeom_log.h"
+/* GeographicLib */
+#include "geodesic.h"
+
/**
* Initialize spheroid object based on major and minor axis
*/
s->radius = (2.0 * a + b ) / 3.0;
}
+#ifdef USE_PRE22GEODESIC
static double spheroid_mu2(double alpha, const SPHEROID *s)
{
double b2 = POW2(s->b);
{
return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
}
+#endif /* def USE_PRE22GEODESIC */
+
+
+#ifndef USE_PRE22GEODESIC
+
+/**
+* Computes the shortest distance along the surface of the spheroid
+* between two points, using the inverse geodesic problem from
+* GeographicLib (Karney 2013).
+*
+* @param a - location of first point
+* @param b - location of second point
+* @param s - spheroid to calculate on
+* @return spheroidal distance between a and b in spheroid units
+*/
+double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
+{
+ struct geod_geodesic gd;
+ geod_init(&gd, spheroid->a, spheroid->f);
+ double lat1 = a->lat * 180.0 / M_PI;
+ double lon1 = a->lon * 180.0 / M_PI;
+ double lat2 = b->lat * 180.0 / M_PI;
+ double lon2 = b->lon * 180.0 / M_PI;
+ double s12; /* return distance */
+ geod_inverse(&gd, lat1, lon1, lat2, lon2, &s12, 0, 0);
+ return s12;
+}
+
+/**
+* Computes the forward azimuth of the geodesic joining two points on
+* the spheroid, using the inverse geodesic problem (Karney 2013).
+*
+* @param r - location of first point
+* @param s - location of second point
+* @return azimuth of line joining r to s (but not reverse)
+*/
+double spheroid_direction(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
+{
+ struct geod_geodesic gd;
+ geod_init(&gd, spheroid->a, spheroid->f);
+ double lat1 = a->lat * 180.0 / M_PI;
+ double lon1 = a->lon * 180.0 / M_PI;
+ double lat2 = b->lat * 180.0 / M_PI;
+ double lon2 = b->lon * 180.0 / M_PI;
+ double azi1; /* return azimuth */
+ geod_inverse(&gd, lat1, lon1, lat2, lon2, 0, &azi1, 0);
+ return azi1 * M_PI / 180.0;
+}
+
+/**
+* Given a location, an azimuth and a distance, computes the location of
+* the projected point. Using the direct geodesic problem from
+* GeographicLib (Karney 2013).
+*
+* @param r - location of first point
+* @param distance - distance in meters
+* @param azimuth - azimuth in radians
+* @return g - location of projected point
+*/
+int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
+{
+ struct geod_geodesic gd;
+ geod_init(&gd, spheroid->a, spheroid->f);
+ double lat1 = r->lat * 180.0 / M_PI;
+ double lon1 = r->lon * 180.0 / M_PI;
+ double lat2, lon2; /* return projected position */
+ geod_direct(&gd, lat1, lon1, azimuth * 180.0 / M_PI, distance, &lat2, &lon2, 0);
+ g->lat = lat2 * M_PI / 180.0;
+ g->lon = lon2 * M_PI / 180.0;
+ return LW_SUCCESS;
+}
+
+
+static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
+{
+ /* Return zero on non-sensical inputs */
+ if ( ! pa || pa->npoints < 4 )
+ return 0.0;
+
+ struct geod_geodesic gd;
+ geod_init(&gd, spheroid->a, spheroid->f);
+ struct geod_polygon poly;
+ geod_polygon_init(&poly, 0);
+ int i;
+ double area; /* returned polygon area */
+ POINT2D p; /* long/lat units are degrees */
+
+ /* Pass points from point array; don't close the linearring */
+ for ( i = 0; i < pa->npoints - 1; i++ )
+ {
+ getPoint2d_p(pa, i, &p);
+ geod_polygon_addpoint(&gd, &poly, p.y, p.x);
+ LWDEBUGF(4, "geod_polygon_addpoint %d: %.12g %.12g", i, p.y, p.x);
+ }
+ i = geod_polygon_compute(&gd, &poly, 0, 1, &area, 0);
+ if ( i != pa->npoints - 1 )
+ {
+ lwerror("ptarray_area_spheroid: different number of points %d vs %d",
+ i, pa->npoints - 1);
+ }
+ LWDEBUGF(4, "geod_polygon_compute area: %.12g", area);
+ return fabs(area);
+}
+
+/* Above use GeographicLib */
+#else /* ndef USE_PRE22GEODESIC */
+/* Below use pre-version 2.2 geodesic functions */
/**
* Computes the shortest distance along the surface of the spheroid
{
azimuth = azimuth + M_PI * 2.0;
}
- if (azimuth > (PI * 2.0))
+ if (azimuth > (M_PI * 2.0))
{
azimuth = azimuth - M_PI * 2.0;
}
}
return fabs(area);
}
+#endif /* else USE_PRE22GEODESIC */
+
/**
* Calculate the area of an LWGEOM. Anything except POLYGON, MULTIPOLYGON
* and GEOMETRYCOLLECTION return zero immediately. Multi's recurse, polygons
#include "geography_measurement_trees.h" /* For circ_tree caching */
#include "lwgeom_transform.h" /* For SRID functions */
+#ifdef USE_PRE22GEODESIC
+/* round to 100 nm precision */
+#define INVMINDIST 1.0e9
+#else
+/* round to 10 nm precision */
+#define INVMINDIST 1.0e8
+#endif
+
Datum geography_distance(PG_FUNCTION_ARGS);
Datum geography_distance_uncached(PG_FUNCTION_ARGS);
Datum geography_distance_tree(PG_FUNCTION_ARGS);
PG_FREE_IF_COPY(g1, 0);
PG_FREE_IF_COPY(g2, 1);
+ /* Knock off any funny business at the nanometer level, ticket #2168 */
+ distance = round(distance * INVMINDIST) / INVMINDIST;
+
/* Something went wrong, negative return... should already be eloged, return NULL */
if ( distance < 0.0 )
{
PG_RETURN_NULL();
}
- /* Knock off any funny business at the micrometer level, ticket #2168 */
- distance = round(distance * 10e8) / 10e8;
-
PG_RETURN_FLOAT8(distance);
}
else
lwgeom_calculate_gbox_geodetic(lwgeom, &gbox);
+#ifdef USE_PRE22GEODESIC
/* Test for cases that are currently not handled by spheroid code */
if ( use_spheroid )
{
if ( gbox.zmax > 0.0 && gbox.zmin < 0.0 )
use_spheroid = LW_FALSE;
}
+#endif /* ifdef USE_PRE22GEODESIC */
/* User requests spherical calculation, turn our spheroid into a sphere */
if ( ! use_spheroid )
return sphere->b * (A * (sigma - dsigma));
}
+
/*
* Find the "length of a geometry"
* length2d_spheroid(point, sphere) = 0
* length2d_spheroid(line, sphere) = length of line
* length2d_spheroid(polygon, sphere) = 0
- * -- could make sense to return sum(ring perimeter)
+ * -- could make sense to return sum(ring perimeter)
* uses ellipsoidal math to find the distance
* x's are longitude, and y's are latitude - both in decimal degrees
*/
PG_FUNCTION_INFO_V1(LWGEOM_length2d_ellipsoid);
Datum LWGEOM_length2d_ellipsoid(PG_FUNCTION_ARGS)
{
- Datum geom2d = DirectFunctionCall1(LWGEOM_force_2d, PG_GETARG_DATUM(0));
- PG_RETURN_DATUM(DirectFunctionCall2(LWGEOM_length_ellipsoid_linestring,
- geom2d, PG_GETARG_DATUM(1)));
+ GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
+ SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
+ LWGEOM *lwgeom = lwgeom_from_gserialized(geom);
+ double dist = lwgeom_length_spheroid(lwgeom, sphere);
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 0);
+ PG_RETURN_FLOAT8(dist);
}
R = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * sin2_lat));
/* 90 - lat1, but in radians */
- S = R * sin(M_PI/2.0-lat1) ;
+ S = R * sin(M_PI_2 - lat1) ;
deltaX = long2 - long1; /* in rads */
deltaY = lat2 - lat1; /* in rads */
curver_text text;
BEGIN
--
- -- Raises a DEBUG if it was deprecated in this version,
+ -- Raises a NOTICE if it was deprecated in this version,
-- a WARNING if in a previous version (only up to minor version checked)
--
curver_text := POSTGIS_LIB_VERSION;
emptyCollectionArea|0
spheroidLength1_deprecated|85204.52077
spheroidLength1|85204.52077
-length2d_spheroid|0
+length2d_spheroid|100
length_spheroid|100
projects / postgis / commitdiff