PostGIS  2.5.2dev-r@@SVN_REVISION@@

◆ geography_centroid_from_mpoly()

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 311 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().

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
double x
Definition: liblwgeom.h:348
uint32_t ngeoms
Definition: liblwgeom.h:496
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
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
uint32_t nrings
Definition: liblwgeom.h:457
unsigned int uint32_t
Definition: uthash.h:78
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:63
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
POINTARRAY ** rings
Definition: liblwgeom.h:459
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
LWPOLY * lwpoly_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoly.c:161
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:354
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...
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
uint32_t npoints
Definition: liblwgeom.h:373
Here is the call graph for this function:
Here is the caller graph for this function: