From 9726c47706662604f51411f36fd521d7f0af5de9 Mon Sep 17 00:00:00 2001 From: Regina Obe Date: Fri, 4 Aug 2017 16:23:31 +0000 Subject: [PATCH] missed file commit for ST_Centroid Geography support references #2951 git-svn-id: http://svn.osgeo.org/postgis/trunk@15519 b70326c6-7e19-0410-871a-916f4a2858ee --- postgis/geography_centroid.c | 386 +++++++++++++++++++++++++++++++++++ 1 file changed, 386 insertions(+) create mode 100644 postgis/geography_centroid.c diff --git a/postgis/geography_centroid.c b/postgis/geography_centroid.c new file mode 100644 index 000000000..c7710c84f --- /dev/null +++ b/postgis/geography_centroid.c @@ -0,0 +1,386 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * PostGIS is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * PostGIS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PostGIS. If not, see . + * + ********************************************************************** + * + * Copyright (C) 2017 Danny Götte + * + **********************************************************************/ + +#include "postgres.h" + +#include "../postgis_config.h" + +#include + +#include "liblwgeom.h" /* For standard geometry types. */ +#include "lwgeom_pg.h" /* For pg macros. */ +#include "lwgeom_transform.h" /* For SRID functions */ + +Datum geography_centroid(PG_FUNCTION_ARGS); + +/* internal functions */ +LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size); +LWPOINT* geography_centroid_from_mline(const LWMLINE* mline, SPHEROID* s); +LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s); +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); +POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat); + +/** + * geography_centroid(GSERIALIZED *g) + * returns centroid as point + */ +PG_FUNCTION_INFO_V1(geography_centroid); +Datum geography_centroid(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom = NULL; + LWGEOM *lwgeom_out = NULL; + LWPOINT *lwpoint_out = NULL; + GSERIALIZED *g = NULL; + GSERIALIZED *g_out = NULL; + uint32_t srid; + bool use_spheroid = true; + SPHEROID s; + uint32_t type; + + /* Get our geometry object loaded into memory. */ + g = PG_GETARG_GSERIALIZED_P(0); + lwgeom = lwgeom_from_gserialized(g); + + if (g == NULL) + { + PG_RETURN_NULL(); + } + + srid = lwgeom_get_srid(lwgeom); + + /* on empty input, return empty output */ + if (gserialized_is_empty(g)) + { + LWCOLLECTION* empty = lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0); + lwgeom_out = lwcollection_as_lwgeom(empty); + lwgeom_set_geodetic(lwgeom_out, true); + g_out = gserialized_from_lwgeom(lwgeom_out, 0); + PG_RETURN_POINTER(g_out); + } + + /* Initialize spheroid */ + spheroid_init_from_srid(fcinfo, srid, &s); + + /* Set to sphere if requested */ + use_spheroid = PG_GETARG_BOOL(1); + if ( ! use_spheroid ) + s.a = s.b = s.radius; + + type = gserialized_get_type(g); + + switch (type) + { + + case POINTTYPE: + { + /* centroid of a point is itself */ + PG_RETURN_POINTER(g); + break; + } + + case MULTIPOINTTYPE: + { + LWMPOINT* mpoints = lwgeom_as_lwmpoint(lwgeom); + + /* average between all points */ + uint32_t size = mpoints->ngeoms; + POINT3DM points[size]; + + uint32_t i; + for (i = 0; i < size; i++) { + points[i].x = lwpoint_get_x(mpoints->geoms[i]); + points[i].y = lwpoint_get_y(mpoints->geoms[i]); + points[i].m = 1; + } + + lwpoint_out = geography_centroid_from_wpoints(srid, points, size); + break; + } + + case LINETYPE: + { + LWLINE* line = lwgeom_as_lwline(lwgeom); + + /* reuse mline function */ + LWMLINE* mline = lwmline_construct_empty(srid, 0, 0); + lwmline_add_lwline(mline, line); + + lwpoint_out = geography_centroid_from_mline(mline, &s); + lwmline_free(mline); + break; + } + + case MULTILINETYPE: + { + LWMLINE* mline = lwgeom_as_lwmline(lwgeom); + lwpoint_out = geography_centroid_from_mline(mline, &s); + break; + } + + case POLYGONTYPE: + { + LWPOLY* poly = lwgeom_as_lwpoly(lwgeom); + + /* reuse mpoly function */ + LWMPOLY* mpoly = lwmpoly_construct_empty(srid, 0, 0); + lwmpoly_add_lwpoly(mpoly, poly); + + lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s); + lwmpoly_free(mpoly); + break; + } + + case MULTIPOLYGONTYPE: + { + LWMPOLY* mpoly = lwgeom_as_lwmpoly(lwgeom); + lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s); + break; + } + default: + elog(ERROR, "ST_Centroid(geography) unhandled geography type"); + PG_RETURN_NULL(); + } + + PG_FREE_IF_COPY(g, 0); + + lwgeom_out = lwpoint_as_lwgeom(lwpoint_out); + lwgeom_set_geodetic(lwgeom_out, true); + g_out = gserialized_from_lwgeom(lwgeom_out, 0); + + PG_RETURN_POINTER(g_out); +} + + +/** + * Convert lat-lon-points to x-y-z-coordinates, calculate a weighted average + * point and return lat-lon-coordinated + */ +LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size) +{ + double_t x_sum = 0; + double_t y_sum = 0; + double_t z_sum = 0; + double_t weight_sum = 0; + + double_t weight = 1; + POINT3D* point; + + uint32_t i; + for (i = 0; i < size; i++ ) + { + point = lonlat_to_cart(points[i].x, points[i].y); + weight = points[i].m; + + x_sum += point->x * weight; + y_sum += point->y * weight; + z_sum += point->z * weight; + + weight_sum += weight; + + lwfree(point); + } + + return cart_to_lwpoint(x_sum, y_sum, z_sum, weight_sum, srid); +} + +POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat) +{ + double_t lat, lon; + double_t sin_lat; + + POINT3D* point = lwalloc(sizeof(POINT3D));; + + // prepare coordinate for trigonometric functions from [-90, 90] -> [0, pi] + lat = (raw_lat + 90) / 180 * M_PI; + + // prepare coordinate for trigonometric functions from [-180, 180] -> [-pi, pi] + lon = raw_lon / 180 * M_PI; + + /* calculate value only once */ + sin_lat = sinl(lat); + + /* convert to 3D cartesian coordinates */ + point->x = sin_lat * cosl(lon); + point->y = sin_lat * sinl(lon); + point->z = cosl(lat); + + return point; +} + +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) +{ + double_t x = x_sum / weight_sum; + double_t y = y_sum / weight_sum; + double_t z = z_sum / weight_sum; + + /* x-y-z vector length */ + double_t r = sqrtl(powl(x, 2) + powl(y, 2) + powl(z, 2)); + + double_t lon = atan2l(y, x) * 180 / M_PI; + double_t lat = acosl(z / r) * 180 / M_PI - 90; + + return lwpoint_make2d(srid, lon, lat); +} + +/** + * Split lines into segments and calculate with middle of segment as weighted + * point. + */ +LWPOINT* geography_centroid_from_mline(const LWMLINE* mline, SPHEROID* s) +{ + double_t tolerance = 0.0; + uint32_t size = 0; + uint32_t i, k; + + /* get total number of points */ + for (i = 0; i < mline->ngeoms; i++) { + size += (mline->geoms[i]->points->npoints - 1) * 2; + } + + POINT3DM points[size]; + uint32_t j = 0; + + for (i = 0; i < mline->ngeoms; i++) { + LWLINE* line = mline->geoms[i]; + + /* add both points of line segment as weighted point */ + for (k = 0; k < line->points->npoints - 1; k++) { + const POINT2D* p1 = getPoint2d_cp(line->points, k); + const POINT2D* p2 = getPoint2d_cp(line->points, k+1); + + /* use line-segment length as weight */ + LWPOINT* lwp1 = lwpoint_make2d(mline->srid, p1->x, p1->y); + LWPOINT* lwp2 = lwpoint_make2d(mline->srid, p2->x, p2->y); + LWGEOM* lwgeom1 = lwpoint_as_lwgeom(lwp1); + LWGEOM* lwgeom2 = lwpoint_as_lwgeom(lwp2); + lwgeom_set_geodetic(lwgeom1, LW_TRUE); + lwgeom_set_geodetic(lwgeom2, LW_TRUE); + + /* use point distance as weight */ + double_t weight = lwgeom_distance_spheroid(lwgeom1, lwgeom2, s, tolerance); + + points[j].x = p1->x; + points[j].y = p1->y; + points[j].m = weight; + j++; + + points[j].x = p2->x; + points[j].y = p2->y; + points[j].m = weight; + j++; + + lwgeom_free(lwgeom1); + lwgeom_free(lwgeom2); + } + } + + return geography_centroid_from_wpoints(mline->srid, points, size); +} + + +/** + * Split polygons into triangles and use centroid of the triangle with the + * triangle area as weight to calculate the centroid of a (multi)polygon. + */ +LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s) +{ + uint32_t size = 0; + uint32_t i, ir, ip; + for (ip = 0; ip < mpoly->ngeoms; ip++) { + for (ir = 0; ir < mpoly->geoms[ip]->nrings; ir++) { + size += mpoly->geoms[ip]->rings[ir]->npoints; + } + } + + POINT3DM points[size]; + uint32_t j = 0; + + GBOX gbox; + + /* use first point as reference to create triangles */ + const POINT4D* reference_point = (const POINT4D*) getPoint2d_cp(mpoly->geoms[0]->rings[0], 0); + + for (ip = 0; ip < mpoly->ngeoms; ip++) { + LWPOLY* poly = mpoly->geoms[ip]; + + for (ir = 0; ir < poly->nrings; ir++) { + POINTARRAY* ring = poly->rings[ir]; + + /* split into triangles (two points + reference point) */ + for (i = 0; i < ring->npoints - 1; i++) { + const POINT4D* p1 = (const POINT4D*) getPoint2d_cp(ring, i); + const POINT4D* p2 = (const POINT4D*) getPoint2d_cp(ring, i+1); + + POINTARRAY* pa = ptarray_construct_empty(0, 0, 4); + ptarray_insert_point(pa, p1, 0); + ptarray_insert_point(pa, p2, 1); + ptarray_insert_point(pa, reference_point, 2); + ptarray_insert_point(pa, p1, 3); + + LWPOLY* poly_tri = lwpoly_construct_empty(mpoly->srid, 0, 0); + lwpoly_add_ring(poly_tri, pa); + + LWGEOM* geom_tri = lwpoly_as_lwgeom(poly_tri); + lwgeom_set_geodetic(geom_tri, LW_TRUE); + + /* Calculate the weight of the triangle. If counter clockwise, + * the weight is negative (e.g. for holes in polygons) + */ + + double_t weight; + if ( use_spheroid ) + weight = lwgeom_area_spheroid(geom_tri, s); + else + weight = lwgeom_area_sphere(geom_tri, s); + + + POINT3DM triangle[3]; + triangle[0].x = p1->x; + triangle[0].y = p1->y; + triangle[0].m = 1; + + triangle[1].x = p2->x; + triangle[1].y = p2->y; + triangle[1].m = 1; + + triangle[2].x = reference_point->x; + triangle[2].y = reference_point->y; + triangle[2].m = 1; + + /* get center of triangle */ + LWPOINT* tri_centroid = geography_centroid_from_wpoints(mpoly->srid, triangle, 3); + + points[j].x = lwpoint_get_x(tri_centroid); + points[j].y = lwpoint_get_y(tri_centroid); + points[j].m = weight; + j++; + + lwpoint_free(tri_centroid); + lwgeom_free(geom_tri); + } + } + } + + return geography_centroid_from_wpoints(mpoly->srid, points, size); +} -- 2.40.0