PostGIS  2.5.2dev-r@@SVN_REVISION@@
geography_centroid.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2017 Danny Götte <danny.goette@fem.tu-ilmenau.de>
22  *
23  **********************************************************************/
24 
25 #include "postgres.h"
26 
27 #include "../postgis_config.h"
28 
29 #include <math.h>
30 
31 #include "liblwgeom.h" /* For standard geometry types. */
32 #include "lwgeom_pg.h" /* For pg macros. */
33 #include "lwgeom_transform.h" /* For SRID functions */
34 
35 Datum geography_centroid(PG_FUNCTION_ARGS);
36 
37 /* internal functions */
38 LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size);
40 LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s);
41 LWPOINT* cart_to_lwpoint(const double_t x_sum, const double_t y_sum, const double_t z_sum, const double_t weight_sum, const uint32_t srid);
42 POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat);
43 
49 Datum geography_centroid(PG_FUNCTION_ARGS)
50 {
51  LWGEOM *lwgeom = NULL;
52  LWGEOM *lwgeom_out = NULL;
53  LWPOINT *lwpoint_out = NULL;
54  GSERIALIZED *g = NULL;
55  GSERIALIZED *g_out = NULL;
56  uint32_t srid;
57  bool use_spheroid = true;
58  SPHEROID s;
59  uint32_t type;
60 
61  /* Get our geometry object loaded into memory. */
62  g = PG_GETARG_GSERIALIZED_P(0);
63  lwgeom = lwgeom_from_gserialized(g);
64 
65  if (g == NULL)
66  {
67  PG_RETURN_NULL();
68  }
69 
70  srid = lwgeom_get_srid(lwgeom);
71 
72  /* on empty input, return empty output */
73  if (gserialized_is_empty(g))
74  {
76  lwgeom_out = lwcollection_as_lwgeom(empty);
77  lwgeom_set_geodetic(lwgeom_out, true);
78  g_out = gserialized_from_lwgeom(lwgeom_out, 0);
79  PG_RETURN_POINTER(g_out);
80  }
81 
82  /* Initialize spheroid */
83  spheroid_init_from_srid(fcinfo, srid, &s);
84 
85  /* Set to sphere if requested */
86  use_spheroid = PG_GETARG_BOOL(1);
87  if ( ! use_spheroid )
88  s.a = s.b = s.radius;
89 
90  type = gserialized_get_type(g);
91 
92  switch (type)
93  {
94 
95  case POINTTYPE:
96  {
97  /* centroid of a point is itself */
98  PG_RETURN_POINTER(g);
99  break;
100  }
101 
102  case MULTIPOINTTYPE:
103  {
104  LWMPOINT* mpoints = lwgeom_as_lwmpoint(lwgeom);
105 
106  /* average between all points */
107  uint32_t size = mpoints->ngeoms;
108  POINT3DM* points = palloc(size*sizeof(POINT3DM));
109 
110  uint32_t i;
111  for (i = 0; i < size; i++) {
112  points[i].x = lwpoint_get_x(mpoints->geoms[i]);
113  points[i].y = lwpoint_get_y(mpoints->geoms[i]);
114  points[i].m = 1;
115  }
116 
117  lwpoint_out = geography_centroid_from_wpoints(srid, points, size);
118  pfree(points);
119  break;
120  }
121 
122  case LINETYPE:
123  {
124  LWLINE* line = lwgeom_as_lwline(lwgeom);
125 
126  /* reuse mline function */
127  LWMLINE* mline = lwmline_construct_empty(srid, 0, 0);
128  lwmline_add_lwline(mline, line);
129 
130  lwpoint_out = geography_centroid_from_mline(mline, &s);
131  lwmline_free(mline);
132  break;
133  }
134 
135  case MULTILINETYPE:
136  {
137  LWMLINE* mline = lwgeom_as_lwmline(lwgeom);
138  lwpoint_out = geography_centroid_from_mline(mline, &s);
139  break;
140  }
141 
142  case POLYGONTYPE:
143  {
144  LWPOLY* poly = lwgeom_as_lwpoly(lwgeom);
145 
146  /* reuse mpoly function */
147  LWMPOLY* mpoly = lwmpoly_construct_empty(srid, 0, 0);
148  lwmpoly_add_lwpoly(mpoly, poly);
149 
150  lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
151  lwmpoly_free(mpoly);
152  break;
153  }
154 
155  case MULTIPOLYGONTYPE:
156  {
157  LWMPOLY* mpoly = lwgeom_as_lwmpoly(lwgeom);
158  lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
159  break;
160  }
161  default:
162  elog(ERROR, "ST_Centroid(geography) unhandled geography type");
163  PG_RETURN_NULL();
164  }
165 
166  PG_FREE_IF_COPY(g, 0);
167 
168  lwgeom_out = lwpoint_as_lwgeom(lwpoint_out);
169  lwgeom_set_geodetic(lwgeom_out, true);
170  g_out = gserialized_from_lwgeom(lwgeom_out, 0);
171 
172  PG_RETURN_POINTER(g_out);
173 }
174 
175 
180 LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size)
181 {
182  double_t x_sum = 0;
183  double_t y_sum = 0;
184  double_t z_sum = 0;
185  double_t weight_sum = 0;
186 
187  double_t weight = 1;
188  POINT3D* point;
189 
190  uint32_t i;
191  for (i = 0; i < size; i++ )
192  {
193  point = lonlat_to_cart(points[i].x, points[i].y);
194  weight = points[i].m;
195 
196  x_sum += point->x * weight;
197  y_sum += point->y * weight;
198  z_sum += point->z * weight;
199 
200  weight_sum += weight;
201 
202  lwfree(point);
203  }
204 
205  return cart_to_lwpoint(x_sum, y_sum, z_sum, weight_sum, srid);
206 }
207 
208 POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat)
209 {
210  double_t lat, lon;
211  double_t sin_lat;
212 
213  POINT3D* point = lwalloc(sizeof(POINT3D));;
214 
215  // prepare coordinate for trigonometric functions from [-90, 90] -> [0, pi]
216  lat = (raw_lat + 90) / 180 * M_PI;
217 
218  // prepare coordinate for trigonometric functions from [-180, 180] -> [-pi, pi]
219  lon = raw_lon / 180 * M_PI;
220 
221  /* calculate value only once */
222  sin_lat = sinl(lat);
223 
224  /* convert to 3D cartesian coordinates */
225  point->x = sin_lat * cosl(lon);
226  point->y = sin_lat * sinl(lon);
227  point->z = cosl(lat);
228 
229  return point;
230 }
231 
232 LWPOINT* cart_to_lwpoint(const double_t x_sum, const double_t y_sum, const double_t z_sum, const double_t weight_sum, const uint32_t srid)
233 {
234  double_t x = x_sum / weight_sum;
235  double_t y = y_sum / weight_sum;
236  double_t z = z_sum / weight_sum;
237 
238  /* x-y-z vector length */
239  double_t r = sqrtl(powl(x, 2) + powl(y, 2) + powl(z, 2));
240 
241  double_t lon = atan2l(y, x) * 180 / M_PI;
242  double_t lat = acosl(z / r) * 180 / M_PI - 90;
243 
244  return lwpoint_make2d(srid, lon, lat);
245 }
246 
252 {
253  double_t tolerance = 0.0;
254  uint32_t size = 0;
255  uint32_t i, k, j = 0;
256  POINT3DM* points;
257  LWPOINT* result;
258 
259  /* get total number of points */
260  for (i = 0; i < mline->ngeoms; i++) {
261  size += (mline->geoms[i]->points->npoints - 1) * 2;
262  }
263 
264  points = palloc(size*sizeof(POINT3DM));
265 
266  for (i = 0; i < mline->ngeoms; i++) {
267  LWLINE* line = mline->geoms[i];
268 
269  /* add both points of line segment as weighted point */
270  for (k = 0; k < line->points->npoints - 1; k++) {
271  const POINT2D* p1 = getPoint2d_cp(line->points, k);
272  const POINT2D* p2 = getPoint2d_cp(line->points, k+1);
273  double_t weight;
274 
275  /* use line-segment length as weight */
276  LWPOINT* lwp1 = lwpoint_make2d(mline->srid, p1->x, p1->y);
277  LWPOINT* lwp2 = lwpoint_make2d(mline->srid, p2->x, p2->y);
278  LWGEOM* lwgeom1 = lwpoint_as_lwgeom(lwp1);
279  LWGEOM* lwgeom2 = lwpoint_as_lwgeom(lwp2);
280  lwgeom_set_geodetic(lwgeom1, LW_TRUE);
281  lwgeom_set_geodetic(lwgeom2, LW_TRUE);
282 
283  /* use point distance as weight */
284  weight = lwgeom_distance_spheroid(lwgeom1, lwgeom2, s, tolerance);
285 
286  points[j].x = p1->x;
287  points[j].y = p1->y;
288  points[j].m = weight;
289  j++;
290 
291  points[j].x = p2->x;
292  points[j].y = p2->y;
293  points[j].m = weight;
294  j++;
295 
296  lwgeom_free(lwgeom1);
297  lwgeom_free(lwgeom2);
298  }
299  }
300 
301  result = geography_centroid_from_wpoints(mline->srid, points, size);
302  pfree(points);
303  return result;
304 }
305 
306 
311 LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s)
312 {
313  uint32_t size = 0;
314  uint32_t i, ir, ip, j = 0;
315  POINT3DM* points;
316  POINT4D* reference_point = NULL;
317  LWPOINT* result = NULL;
318 
319  for (ip = 0; ip < mpoly->ngeoms; ip++) {
320  for (ir = 0; ir < mpoly->geoms[ip]->nrings; ir++) {
321  size += mpoly->geoms[ip]->rings[ir]->npoints - 1;
322  }
323  }
324 
325  points = palloc(size*sizeof(POINT3DM));
326 
327 
328  /* use first point as reference to create triangles */
329  reference_point = (POINT4D*) getPoint2d_cp(mpoly->geoms[0]->rings[0], 0);
330 
331  for (ip = 0; ip < mpoly->ngeoms; ip++) {
332  LWPOLY* poly = mpoly->geoms[ip];
333 
334  for (ir = 0; ir < poly->nrings; ir++) {
335  POINTARRAY* ring = poly->rings[ir];
336 
337  /* split into triangles (two points + reference point) */
338  for (i = 0; i < ring->npoints - 1; i++) {
339  const POINT4D* p1 = (const POINT4D*) getPoint2d_cp(ring, i);
340  const POINT4D* p2 = (const POINT4D*) getPoint2d_cp(ring, i+1);
341  LWPOLY* poly_tri;
342  LWGEOM* geom_tri;
343  double_t weight;
344  POINT3DM triangle[3];
345  LWPOINT* tri_centroid;
346 
347  POINTARRAY* pa = ptarray_construct_empty(0, 0, 4);
348  ptarray_insert_point(pa, p1, 0);
349  ptarray_insert_point(pa, p2, 1);
350  ptarray_insert_point(pa, reference_point, 2);
351  ptarray_insert_point(pa, p1, 3);
352 
353  poly_tri = lwpoly_construct_empty(mpoly->srid, 0, 0);
354  lwpoly_add_ring(poly_tri, pa);
355 
356  geom_tri = lwpoly_as_lwgeom(poly_tri);
357  lwgeom_set_geodetic(geom_tri, LW_TRUE);
358 
359  /* Calculate the weight of the triangle. If counter clockwise,
360  * the weight is negative (e.g. for holes in polygons)
361  */
362 
363  if ( use_spheroid )
364  weight = lwgeom_area_spheroid(geom_tri, s);
365  else
366  weight = lwgeom_area_sphere(geom_tri, s);
367 
368 
369  triangle[0].x = p1->x;
370  triangle[0].y = p1->y;
371  triangle[0].m = 1;
372 
373  triangle[1].x = p2->x;
374  triangle[1].y = p2->y;
375  triangle[1].m = 1;
376 
377  triangle[2].x = reference_point->x;
378  triangle[2].y = reference_point->y;
379  triangle[2].m = 1;
380 
381  /* get center of triangle */
382  tri_centroid = geography_centroid_from_wpoints(mpoly->srid, triangle, 3);
383 
384  points[j].x = lwpoint_get_x(tri_centroid);
385  points[j].y = lwpoint_get_y(tri_centroid);
386  points[j].m = weight;
387  j++;
388 
389  lwpoint_free(tri_centroid);
390  lwgeom_free(geom_tri);
391  }
392  }
393  }
394  result = geography_centroid_from_wpoints(mpoly->srid, points, size);
395  pfree(points);
396  return result;
397 }
double x
Definition: liblwgeom.h:354
#define LINETYPE
Definition: liblwgeom.h:85
uint32_t gserialized_get_type(const GSERIALIZED *s)
Extract the geometry type from the serialized form (it hides in the anonymous data area...
Definition: g_serialized.c:86
uint32_t ngeoms
Definition: liblwgeom.h:470
LWPOINT * cart_to_lwpoint(const double_t x_sum, const double_t y_sum, const double_t z_sum, const double_t weight_sum, const uint32_t srid)
char * r
Definition: cu_in_wkt.c:24
void lwmline_free(LWMLINE *mline)
Definition: lwmline.c:112
void lwfree(void *mem)
Definition: lwutil.c:244
double x
Definition: liblwgeom.h:348
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
double y
Definition: liblwgeom.h:342
uint32_t ngeoms
Definition: liblwgeom.h:496
#define POLYGONTYPE
Definition: liblwgeom.h:86
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
double b
Definition: liblwgeom.h:316
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
Datum geography_centroid(PG_FUNCTION_ARGS)
LWMPOLY * lwmpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoly.c:40
double radius
Definition: liblwgeom.h:320
POINT3D * lonlat_to_cart(const double_t raw_lon, const double_t raw_lat)
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:916
double x
Definition: liblwgeom.h:342
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:206
LWMLINE * lwmline_add_lwline(LWMLINE *mobj, const LWLINE *obj)
Definition: lwmline.c:46
LWMLINE * lwmline_construct_empty(int srid, char hasz, char hasm)
Definition: lwmline.c:38
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
uint32_t nrings
Definition: liblwgeom.h:457
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: g_serialized.c:179
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:251
double z
Definition: liblwgeom.h:342
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:63
void lwmpoly_free(LWMPOLY *mpoly)
Definition: lwmpoly.c:53
double m
Definition: liblwgeom.h:348
LWPOLY ** geoms
Definition: liblwgeom.h:498
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWPOINT * geography_centroid_from_mpoly(const LWMPOLY *mpoly, bool use_spheroid, SPHEROID *s)
Split polygons into triangles and use centroid of the triangle with the triangle area as weight to ca...
LWMLINE * lwgeom_as_lwmline(const LWGEOM *lwgeom)
Definition: lwgeom.c:242
POINTARRAY ** rings
Definition: liblwgeom.h:459
LWPOINT ** geoms
Definition: liblwgeom.h:472
double y
Definition: liblwgeom.h:330
char * s
Definition: cu_in_wkt.c:23
double y
Definition: liblwgeom.h:348
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the sphere.
Definition: lwgeodetic.c:1944
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:952
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
PG_FUNCTION_INFO_V1(geography_centroid)
geography_centroid(GSERIALIZED *g) returns centroid as point
LWPOINT * geography_centroid_from_mline(const LWMLINE *mline, SPHEROID *s)
Split lines into segments and calculate with middle of segment as weighted point. ...
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
uint32_t ngeoms
Definition: liblwgeom.h:483
LWLINE ** geoms
Definition: liblwgeom.h:485
double a
Definition: liblwgeom.h:315
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:161
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:233
int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, uint32_t where)
Insert a point into an existing POINTARRAY.
Definition: ptarray.c:96
type
Definition: ovdump.py:41
double lwpoint_get_y(const LWPOINT *point)
Definition: lwpoint.c:76
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:335
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the spheroid.
Definition: lwspheroid.c:642
int lwpoly_add_ring(LWPOLY *poly, POINTARRAY *pa)
Add a ring, allocating extra space if necessary.
Definition: lwpoly.c:247
double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const SPHEROID *spheroid, double tolerance)
Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid.
Definition: lwgeodetic.c:2100
void * lwalloc(size_t size)
Definition: lwutil.c:229
GSERIALIZED * gserialized_from_lwgeom(LWGEOM *geom, size_t *size)
Allocate a new GSERIALIZED from an LWGEOM.
double y
Definition: liblwgeom.h:354
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
LWMPOLY * lwmpoly_add_lwpoly(LWMPOLY *mobj, const LWPOLY *obj)
Definition: lwmpoly.c:47
int32_t srid
Definition: liblwgeom.h:495
LWPOINT * geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM *points, const uint32_t size)
Convert lat-lon-points to x-y-z-coordinates, calculate a weighted average point and return lat-lon-co...
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
This library is the generic geometry handling section of PostGIS.
int32_t srid
Definition: liblwgeom.h:482
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:374
POINTARRAY * points
Definition: liblwgeom.h:424
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:300
uint32_t npoints
Definition: liblwgeom.h:373