PostGIS  2.5.0dev-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[size];
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  break;
119  }
120 
121  case LINETYPE:
122  {
123  LWLINE* line = lwgeom_as_lwline(lwgeom);
124 
125  /* reuse mline function */
126  LWMLINE* mline = lwmline_construct_empty(srid, 0, 0);
127  lwmline_add_lwline(mline, line);
128 
129  lwpoint_out = geography_centroid_from_mline(mline, &s);
130  lwmline_free(mline);
131  break;
132  }
133 
134  case MULTILINETYPE:
135  {
136  LWMLINE* mline = lwgeom_as_lwmline(lwgeom);
137  lwpoint_out = geography_centroid_from_mline(mline, &s);
138  break;
139  }
140 
141  case POLYGONTYPE:
142  {
143  LWPOLY* poly = lwgeom_as_lwpoly(lwgeom);
144 
145  /* reuse mpoly function */
146  LWMPOLY* mpoly = lwmpoly_construct_empty(srid, 0, 0);
147  lwmpoly_add_lwpoly(mpoly, poly);
148 
149  lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
150  lwmpoly_free(mpoly);
151  break;
152  }
153 
154  case MULTIPOLYGONTYPE:
155  {
156  LWMPOLY* mpoly = lwgeom_as_lwmpoly(lwgeom);
157  lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
158  break;
159  }
160  default:
161  elog(ERROR, "ST_Centroid(geography) unhandled geography type");
162  PG_RETURN_NULL();
163  }
164 
165  PG_FREE_IF_COPY(g, 0);
166 
167  lwgeom_out = lwpoint_as_lwgeom(lwpoint_out);
168  lwgeom_set_geodetic(lwgeom_out, true);
169  g_out = gserialized_from_lwgeom(lwgeom_out, 0);
170 
171  PG_RETURN_POINTER(g_out);
172 }
173 
174 
179 LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size)
180 {
181  double_t x_sum = 0;
182  double_t y_sum = 0;
183  double_t z_sum = 0;
184  double_t weight_sum = 0;
185 
186  double_t weight = 1;
187  POINT3D* point;
188 
189  uint32_t i;
190  for (i = 0; i < size; i++ )
191  {
192  point = lonlat_to_cart(points[i].x, points[i].y);
193  weight = points[i].m;
194 
195  x_sum += point->x * weight;
196  y_sum += point->y * weight;
197  z_sum += point->z * weight;
198 
199  weight_sum += weight;
200 
201  lwfree(point);
202  }
203 
204  return cart_to_lwpoint(x_sum, y_sum, z_sum, weight_sum, srid);
205 }
206 
207 POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat)
208 {
209  double_t lat, lon;
210  double_t sin_lat;
211 
212  POINT3D* point = lwalloc(sizeof(POINT3D));;
213 
214  // prepare coordinate for trigonometric functions from [-90, 90] -> [0, pi]
215  lat = (raw_lat + 90) / 180 * M_PI;
216 
217  // prepare coordinate for trigonometric functions from [-180, 180] -> [-pi, pi]
218  lon = raw_lon / 180 * M_PI;
219 
220  /* calculate value only once */
221  sin_lat = sinl(lat);
222 
223  /* convert to 3D cartesian coordinates */
224  point->x = sin_lat * cosl(lon);
225  point->y = sin_lat * sinl(lon);
226  point->z = cosl(lat);
227 
228  return point;
229 }
230 
231 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)
232 {
233  double_t x = x_sum / weight_sum;
234  double_t y = y_sum / weight_sum;
235  double_t z = z_sum / weight_sum;
236 
237  /* x-y-z vector length */
238  double_t r = sqrtl(powl(x, 2) + powl(y, 2) + powl(z, 2));
239 
240  double_t lon = atan2l(y, x) * 180 / M_PI;
241  double_t lat = acosl(z / r) * 180 / M_PI - 90;
242 
243  return lwpoint_make2d(srid, lon, lat);
244 }
245 
251 {
252  double_t tolerance = 0.0;
253  uint32_t size = 0;
254  uint32_t i, k;
255 
256  /* get total number of points */
257  for (i = 0; i < mline->ngeoms; i++) {
258  size += (mline->geoms[i]->points->npoints - 1) * 2;
259  }
260 
261  POINT3DM points[size];
262  uint32_t j = 0;
263 
264  for (i = 0; i < mline->ngeoms; i++) {
265  LWLINE* line = mline->geoms[i];
266 
267  /* add both points of line segment as weighted point */
268  for (k = 0; k < line->points->npoints - 1; k++) {
269  const POINT2D* p1 = getPoint2d_cp(line->points, k);
270  const POINT2D* p2 = getPoint2d_cp(line->points, k+1);
271 
272  /* use line-segment length as weight */
273  LWPOINT* lwp1 = lwpoint_make2d(mline->srid, p1->x, p1->y);
274  LWPOINT* lwp2 = lwpoint_make2d(mline->srid, p2->x, p2->y);
275  LWGEOM* lwgeom1 = lwpoint_as_lwgeom(lwp1);
276  LWGEOM* lwgeom2 = lwpoint_as_lwgeom(lwp2);
277  lwgeom_set_geodetic(lwgeom1, LW_TRUE);
278  lwgeom_set_geodetic(lwgeom2, LW_TRUE);
279 
280  /* use point distance as weight */
281  double_t weight = lwgeom_distance_spheroid(lwgeom1, lwgeom2, s, tolerance);
282 
283  points[j].x = p1->x;
284  points[j].y = p1->y;
285  points[j].m = weight;
286  j++;
287 
288  points[j].x = p2->x;
289  points[j].y = p2->y;
290  points[j].m = weight;
291  j++;
292 
293  lwgeom_free(lwgeom1);
294  lwgeom_free(lwgeom2);
295  }
296  }
297 
298  return geography_centroid_from_wpoints(mline->srid, points, size);
299 }
300 
301 
306 LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s)
307 {
308  uint32_t size = 0;
309  uint32_t i, ir, ip;
310  for (ip = 0; ip < mpoly->ngeoms; ip++) {
311  for (ir = 0; ir < mpoly->geoms[ip]->nrings; ir++) {
312  size += mpoly->geoms[ip]->rings[ir]->npoints - 1;
313  }
314  }
315 
316  POINT3DM points[size];
317  uint32_t j = 0;
318 
319  /* use first point as reference to create triangles */
320  const POINT4D* reference_point = (const POINT4D*) getPoint2d_cp(mpoly->geoms[0]->rings[0], 0);
321 
322  for (ip = 0; ip < mpoly->ngeoms; ip++) {
323  LWPOLY* poly = mpoly->geoms[ip];
324 
325  for (ir = 0; ir < poly->nrings; ir++) {
326  POINTARRAY* ring = poly->rings[ir];
327 
328  /* split into triangles (two points + reference point) */
329  for (i = 0; i < ring->npoints - 1; i++) {
330  const POINT4D* p1 = (const POINT4D*) getPoint2d_cp(ring, i);
331  const POINT4D* p2 = (const POINT4D*) getPoint2d_cp(ring, i+1);
332 
333  POINTARRAY* pa = ptarray_construct_empty(0, 0, 4);
334  ptarray_insert_point(pa, p1, 0);
335  ptarray_insert_point(pa, p2, 1);
336  ptarray_insert_point(pa, reference_point, 2);
337  ptarray_insert_point(pa, p1, 3);
338 
339  LWPOLY* poly_tri = lwpoly_construct_empty(mpoly->srid, 0, 0);
340  lwpoly_add_ring(poly_tri, pa);
341 
342  LWGEOM* geom_tri = lwpoly_as_lwgeom(poly_tri);
343  lwgeom_set_geodetic(geom_tri, LW_TRUE);
344 
345  /* Calculate the weight of the triangle. If counter clockwise,
346  * the weight is negative (e.g. for holes in polygons)
347  */
348 
349  double_t weight;
350  if ( use_spheroid )
351  weight = lwgeom_area_spheroid(geom_tri, s);
352  else
353  weight = lwgeom_area_sphere(geom_tri, s);
354 
355 
356  POINT3DM triangle[3];
357  triangle[0].x = p1->x;
358  triangle[0].y = p1->y;
359  triangle[0].m = 1;
360 
361  triangle[1].x = p2->x;
362  triangle[1].y = p2->y;
363  triangle[1].m = 1;
364 
365  triangle[2].x = reference_point->x;
366  triangle[2].y = reference_point->y;
367  triangle[2].m = 1;
368 
369  /* get center of triangle */
370  LWPOINT* tri_centroid = geography_centroid_from_wpoints(mpoly->srid, triangle, 3);
371 
372  points[j].x = lwpoint_get_x(tri_centroid);
373  points[j].y = lwpoint_get_y(tri_centroid);
374  points[j].m = weight;
375  j++;
376 
377  lwpoint_free(tri_centroid);
378  lwgeom_free(geom_tri);
379  }
380  }
381  }
382 
383  return geography_centroid_from_wpoints(mpoly->srid, points, size);
384 }
double x
Definition: liblwgeom.h:351
#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:467
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:345
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
double y
Definition: liblwgeom.h:339
uint32_t ngeoms
Definition: liblwgeom.h:493
#define POLYGONTYPE
Definition: liblwgeom.h:86
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
double b
Definition: liblwgeom.h:313
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:1137
#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:317
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:909
double x
Definition: liblwgeom.h:339
LWPOLY * lwgeom_as_lwpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:205
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:319
uint32_t nrings
Definition: liblwgeom.h:454
int gserialized_is_empty(const GSERIALIZED *g)
Check if a GSERIALIZED is empty without deserializing first.
Definition: g_serialized.c:178
LWMPOLY * lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
Definition: lwgeom.c:250
double z
Definition: liblwgeom.h:339
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:327
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:345
LWPOLY ** geoms
Definition: liblwgeom.h:495
#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:241
POINTARRAY ** rings
Definition: liblwgeom.h:456
LWPOINT ** geoms
Definition: liblwgeom.h:469
double y
Definition: liblwgeom.h:327
char * s
Definition: cu_in_wkt.c:23
tuple x
Definition: pixval.py:53
double y
Definition: liblwgeom.h:345
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the geodetic area of a lwgeom on the sphere.
Definition: lwgeodetic.c:1940
void lwgeom_set_geodetic(LWGEOM *geom, int value)
Set the FLAGS geodetic bit on geometry an all sub-geometries and pointlists.
Definition: lwgeom.c:945
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:169
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:480
LWLINE ** geoms
Definition: liblwgeom.h:482
double a
Definition: liblwgeom.h:312
#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:159
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:232
int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, uint32_t where)
Insert a point into an existing POINTARRAY.
Definition: ptarray.c:96
double lwpoint_get_y(const LWPOINT *point)
Definition: lwpoint.c:76
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:334
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:245
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:2096
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:351
#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:492
tuple y
Definition: pixval.py:54
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:479
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:364
POINTARRAY * points
Definition: liblwgeom.h:421
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:299
uint32_t npoints
Definition: liblwgeom.h:370