PostGIS  2.5.0dev-r@@SVN_REVISION@@
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 calculate the centroid of a (multi)polygon.

Definition at line 306 of file geography_centroid.c.

References geography_centroid_from_wpoints(), LWMPOLY::geoms, getPoint2d_cp(), LW_TRUE, lwgeom_area_sphere(), lwgeom_area_spheroid(), lwgeom_free(), lwgeom_set_geodetic(), lwpoint_free(), lwpoint_get_x(), lwpoint_get_y(), lwpoly_add_ring(), lwpoly_as_lwgeom(), lwpoly_construct_empty(), POINT3DM::m, LWMPOLY::ngeoms, POINTARRAY::npoints, LWPOLY::nrings, ptarray_construct_empty(), ptarray_insert_point(), LWPOLY::rings, LWMPOLY::srid, POINT3DM::x, POINT4D::x, POINT3DM::y, and POINT4D::y.

Referenced by geography_centroid().

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
double x
Definition: liblwgeom.h:345
uint32_t ngeoms
Definition: liblwgeom.h:493
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
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:319
uint32_t nrings
Definition: liblwgeom.h:454
unsigned int uint32_t
Definition: uthash.h:78
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:63
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
POINTARRAY ** rings
Definition: liblwgeom.h:456
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
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:159
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
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 y
Definition: liblwgeom.h:351
int32_t srid
Definition: liblwgeom.h:492
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...
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
uint32_t npoints
Definition: liblwgeom.h:370

Here is the call graph for this function:

Here is the caller graph for this function: