38 #define AUTOFIX LW_TRUE
39 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
56 #define GEOS_FREE(...) \
58 geos_destroy((sizeof((void*[]){__VA_ARGS__})/sizeof(void*)), __VA_ARGS__); \
64 lwerror("%s: GEOS Error: %s", __func__, lwgeom_geos_errmsg); \
69 #define GEOS_FAIL_DEBUG() \
72 lwdebug(1, "%s: GEOS Error: %s", __func__, lwgeom_geos_errmsg); \
76 #define GEOS_FREE_AND_FAIL(...) \
78 GEOS_FREE(__VA_ARGS__); \
82 #define GEOS_FREE_AND_FAIL_DEBUG(...) \
85 GEOS_FREE(__VA_ARGS__); \
91 #define RESULT_SRID(...) \
92 (get_result_srid((sizeof((const void*[]){__VA_ARGS__})/sizeof(void*)), __func__, __VA_ARGS__))
101 GEOSGeometry* g = va_arg(ap, GEOSGeometry*);
122 uint32_t size = 0, i;
124 POINT4D point = { 0.0, 0.0, 0.0, 0.0 };
126 LWDEBUG(2,
"ptarray_fromGEOSCoordSeq called");
128 if (!GEOSCoordSeq_getSize(cs, &size))
lwerror(
"Exception thrown");
130 LWDEBUGF(4,
" GEOSCoordSeq size: %d", size);
134 if (!GEOSCoordSeq_getDimensions(cs, &dims))
lwerror(
"Exception thrown");
136 LWDEBUGF(4,
" GEOSCoordSeq dimensions: %d", dims);
139 if (dims > 3) dims = 3;
142 LWDEBUGF(4,
" output dimensions: %d", dims);
146 for (i = 0; i < size; i++)
148 #if POSTGIS_GEOS_VERSION < 38
149 GEOSCoordSeq_getX(cs, i, &(point.
x));
150 GEOSCoordSeq_getY(cs, i, &(point.
y));
151 if (dims >= 3) GEOSCoordSeq_getZ(cs, i, &(point.
z));
154 GEOSCoordSeq_getXYZ(cs, i, &(point.
x), &(point.
y), &(point.
z));
156 GEOSCoordSeq_getXY(cs, i, &(point.
x), &(point.
y));
168 int type = GEOSGeomTypeId(geom);
169 int SRID = GEOSGetSRID(geom);
174 if (want3d && !GEOSHasZ(geom))
176 LWDEBUG(3,
"Geometry has no Z, won't provide one");
182 const GEOSCoordSequence* cs;
184 const GEOSGeometry* g;
189 LWDEBUG(4,
"lwgeom_from_geometry: it's a Point");
190 cs = GEOSGeom_getCoordSeq(geom);
195 case GEOS_LINESTRING:
196 case GEOS_LINEARRING:
197 LWDEBUG(4,
"lwgeom_from_geometry: it's a LineString or LinearRing");
200 cs = GEOSGeom_getCoordSeq(geom);
205 LWDEBUG(4,
"lwgeom_from_geometry: it's a Polygon");
207 ngeoms = GEOSGetNumInteriorRings(geom);
209 g = GEOSGetExteriorRing(geom);
210 cs = GEOSGeom_getCoordSeq(g);
212 for (i = 0; i < ngeoms; i++)
214 g = GEOSGetInteriorRingN(geom, i);
215 cs = GEOSGeom_getCoordSeq(g);
220 case GEOS_MULTIPOINT:
221 case GEOS_MULTILINESTRING:
222 case GEOS_MULTIPOLYGON:
223 case GEOS_GEOMETRYCOLLECTION:
224 LWDEBUG(4,
"lwgeom_from_geometry: it's a Collection or Multi");
226 ngeoms = GEOSGetNumGeometries(geom);
231 for (i = 0; i < ngeoms; i++)
233 g = GEOSGetGeometryN(geom, i);
240 lwerror(
"GEOS2LWGEOM: unknown geometry type: %d",
type);
252 int append_points = 0;
263 lwerror(
"ptarray_to_GEOSCoordSeq called with fix_ring and 0 vertices in ring, cannot fix");
273 if (!(sq = GEOSCoordSeq_create(pa->
npoints + append_points, dims)))
275 lwerror(
"Error creating GEOS Coordinate Sequence");
279 for (i = 0; i < pa->
npoints; i++)
285 LWDEBUGF(4,
"Point: %g,%g,%g", p3d->
x, p3d->
y, p3d->
z);
293 #if POSTGIS_GEOS_VERSION < 38
294 GEOSCoordSeq_setX(sq, i, p2d->
x);
295 GEOSCoordSeq_setY(sq, i, p2d->
y);
296 if (dims == 3) GEOSCoordSeq_setZ(sq, i, p3d->
z);
299 GEOSCoordSeq_setXYZ(sq, i, p2d->
x, p2d->
y, p3d->
z);
301 GEOSCoordSeq_setXY(sq, i, p2d->
x, p2d->
y);
315 for (i = pa->
npoints; i < pa->npoints + append_points; i++)
317 #if POSTGIS_GEOS_VERSION < 38
318 GEOSCoordSeq_setX(sq, i, p2d->
x);
319 GEOSCoordSeq_setY(sq, i, p2d->
y);
321 GEOSCoordSeq_setXY(sq, i, p2d->
x, p2d->
y);
324 if (dims == 3) GEOSCoordSeq_setZ(sq, i, p3d->
z);
331 static inline GEOSGeometry*
337 g = GEOSGeom_createLinearRing(sq);
344 GEOSGeometry* envelope;
346 GEOSCoordSequence* seq = GEOSCoordSeq_create(5, 2);
347 if (!seq)
return NULL;
349 #if POSTGIS_GEOS_VERSION < 38
350 GEOSCoordSeq_setX(seq, 0, box->
xmin);
351 GEOSCoordSeq_setY(seq, 0, box->
ymin);
353 GEOSCoordSeq_setX(seq, 1, box->
xmax);
354 GEOSCoordSeq_setY(seq, 1, box->
ymin);
356 GEOSCoordSeq_setX(seq, 2, box->
xmax);
357 GEOSCoordSeq_setY(seq, 2, box->
ymax);
359 GEOSCoordSeq_setX(seq, 3, box->
xmin);
360 GEOSCoordSeq_setY(seq, 3, box->
ymax);
362 GEOSCoordSeq_setX(seq, 4, box->
xmin);
363 GEOSCoordSeq_setY(seq, 4, box->
ymin);
365 GEOSCoordSeq_setXY(seq, 0, box->
xmin, box->
ymin);
366 GEOSCoordSeq_setXY(seq, 1, box->
xmax, box->
ymin);
367 GEOSCoordSeq_setXY(seq, 2, box->
xmax, box->
ymax);
368 GEOSCoordSeq_setXY(seq, 3, box->
xmin, box->
ymax);
369 GEOSCoordSeq_setXY(seq, 4, box->
xmin, box->
ymin);
372 ring = GEOSGeom_createLinearRing(seq);
375 GEOSCoordSeq_destroy(seq);
379 envelope = GEOSGeom_createPolygon(ring, NULL, 0);
382 GEOSGeom_destroy(ring);
394 GEOSGeom* geoms = NULL;
395 uint32_t ngeoms, i, j;
397 #if LWDEBUG_LEVEL >= 4
413 GEOSGeometry* g =
LWGEOM2GEOS(lwgeom_stroked, autofix);
423 switch (lwgeom->
type)
429 g = GEOSGeom_createEmptyPoint();
432 #if POSTGIS_GEOS_VERSION < 38
434 g = GEOSGeom_createPoint(sq);
439 g = GEOSGeom_createPoint(sq);
444 g = GEOSGeom_createPointFromXY(p->
x, p->
y);
463 g = GEOSGeom_createLineString(sq);
470 g = GEOSGeom_createEmptyPolygon();
474 if (!shell)
return NULL;
475 ngeoms = lwpoly->
nrings - 1;
476 if (ngeoms > 0) geoms =
lwalloc(
sizeof(GEOSGeom) * ngeoms);
478 for (i = 1; i < lwpoly->
nrings; i++)
484 for (k = 0; k < i - 1; k++)
485 GEOSGeom_destroy(geoms[k]);
487 GEOSGeom_destroy(shell);
491 g = GEOSGeom_createPolygon(shell, geoms, ngeoms);
499 g = GEOSGeom_createEmptyPolygon();
506 g = GEOSGeom_createPolygon(shell, NULL, 0);
517 geostype = GEOS_MULTIPOINT;
519 geostype = GEOS_MULTILINESTRING;
521 geostype = GEOS_MULTIPOLYGON;
523 geostype = GEOS_GEOMETRYCOLLECTION;
528 if (ngeoms > 0) geoms =
lwalloc(
sizeof(GEOSGeom) * ngeoms);
531 for (i = 0; i < ngeoms; ++i)
541 for (k = 0; k < j; k++)
542 GEOSGeom_destroy(geoms[k]);
548 g = GEOSGeom_createCollection(geostype, geoms, j);
549 if (ngeoms > 0)
lwfree(geoms);
558 GEOSSetSRID(g, lwgeom->
srid);
560 #if LWDEBUG_LEVEL >= 4
561 wkt = GEOSGeomToWKT(g);
562 LWDEBUGF(4,
"LWGEOM2GEOS: GEOSGeom: %s", wkt);
572 GEOSCoordSequence* seq = GEOSCoordSeq_create(1, 2);
573 GEOSGeometry* geom = NULL;
575 if (!seq)
return NULL;
577 #if POSTGIS_GEOS_VERSION < 38
578 GEOSCoordSeq_setX(seq, 0,
x);
579 GEOSCoordSeq_setY(seq, 0,
y);
581 GEOSCoordSeq_setXY(seq, 0,
x,
y);
584 geom = GEOSGeom_createPoint(seq);
585 if (!geom) GEOSCoordSeq_destroy(seq);
592 GEOSCoordSequence* seq = GEOSCoordSeq_create(2, 2);
593 GEOSGeometry* geom = NULL;
595 if (!seq)
return NULL;
597 #if POSTGIS_GEOS_VERSION < 38
598 GEOSCoordSeq_setX(seq, 0, x1);
599 GEOSCoordSeq_setY(seq, 0, y1);
600 GEOSCoordSeq_setX(seq, 1, x2);
601 GEOSCoordSeq_setY(seq, 1, y2);
603 GEOSCoordSeq_setXY(seq, 0, x1, y1);
604 GEOSCoordSeq_setXY(seq, 1, x2, y2);
607 geom = GEOSGeom_createLineString(seq);
608 if (!geom) GEOSCoordSeq_destroy(seq);
615 const char* ver = GEOSversion();
625 va_start(ap, funcname);
628 for(i = 0; i <
count; i++)
633 lwerror(
"%s: Geometry is null", funcname);
644 lwerror(
"%s: Operation on mixed SRID geometries (%d != %d)", funcname, srid, g->
srid);
667 GEOSSetSRID(g, srid);
705 #if POSTGIS_GEOS_VERSION < 39
706 lwerror(
"Fixed-precision intersection requires GEOS-3.9 or higher");
710 g3 = GEOSIntersectionPrec(g1, g2, prec);
715 g3 = GEOSIntersection(g1, g2);
719 GEOSSetSRID(g3, srid);
745 g3 = GEOSLineMerge(g1);
748 GEOSSetSRID(g3, srid);
783 #if POSTGIS_GEOS_VERSION < 39
784 lwerror(
"Fixed-precision union requires GEOS-3.9 or higher");
788 g3 = GEOSUnaryUnionPrec(g1, prec);
793 g3 = GEOSUnaryUnion(g1);
797 GEOSSetSRID(g3, srid);
819 GEOSGeometry *g1, *g2, *g3;
835 #if POSTGIS_GEOS_VERSION < 39
836 lwerror(
"Fixed-precision difference requires GEOS-3.9 or higher");
840 g3 = GEOSDifferencePrec(g1, g2, prec);
845 g3 = GEOSDifference(g1, g2);
849 GEOSSetSRID(g3, srid);
870 GEOSGeometry *g1, *g2, *g3;
886 #if POSTGIS_GEOS_VERSION < 39
887 lwerror(
"Fixed-precision difference requires GEOS-3.9 or higher");
891 g3 = GEOSSymDifferencePrec(g1, g2, prec);
896 g3 = GEOSSymDifference(g1, g2);
900 GEOSSetSRID(g3, srid);
915 GEOSGeometry *g1, *g3;
929 g3 = GEOSGetCentroid(g1);
932 GEOSSetSRID(g3, srid);
945 #if POSTGIS_GEOS_VERSION < 39
946 lwerror(
"Precision reduction requires GEOS-3.9 or higher");
953 GEOSGeometry *g1, *g3;
964 g3 = GEOSGeom_setPrecision(g1, gridSize, 0);
967 GEOSSetSRID(g3, srid);
984 GEOSGeometry *g1, *g3;
998 g3 = GEOSPointOnSurface(g1);
1001 GEOSSetSRID(g3, srid);
1023 GEOSGeometry *g1, *g2, *g3;
1038 if ( gridSize >= 0) {
1039 #if POSTGIS_GEOS_VERSION < 39
1040 lwerror(
"Fixed-precision union requires GEOS-3.9 or higher");
1044 g3 = GEOSUnionPrec(g1, g2, gridSize);
1049 g3 = GEOSUnion(g1, g2);
1053 GEOSSetSRID(g3, srid);
1066 GEOSGeometry *g1, *g3;
1080 if (!(g3 = GEOSClipByRect(g1, x1, y1, x2, y2)))
1096 #if POSTGIS_GEOS_VERSION < 38
1097 typedef struct Face_t
1099 const GEOSGeometry* geom;
1102 struct Face_t* parent;
1105 static Face* newFace(
const GEOSGeometry* g);
1106 static void delFace(Face* f);
1107 static unsigned int countParens(
const Face* f);
1108 static void findFaceHoles(Face** faces,
int nfaces);
1111 newFace(
const GEOSGeometry* g)
1113 Face* f =
lwalloc(
sizeof(Face));
1115 f->env = GEOSEnvelope(f->geom);
1116 GEOSArea(f->env, &f->envarea);
1122 countParens(
const Face* f)
1124 unsigned int pcount = 0;
1137 GEOSGeom_destroy(f->env);
1142 compare_by_envarea(
const void* g1,
const void* g2)
1144 Face* f1 = *(Face**)g1;
1145 Face* f2 = *(Face**)g2;
1146 double n1 = f1->envarea;
1147 double n2 = f2->envarea;
1149 if (n1 < n2)
return 1;
1150 if (n1 > n2)
return -1;
1156 findFaceHoles(Face** faces,
int nfaces)
1161 qsort(faces, nfaces,
sizeof(Face*), compare_by_envarea);
1162 for (i = 0; i < nfaces; ++i)
1165 int nholes = GEOSGetNumInteriorRings(f->geom);
1166 LWDEBUGF(2,
"Scanning face %d with env area %g and %d holes", i, f->envarea, nholes);
1167 for (h = 0; h < nholes; ++h)
1169 const GEOSGeometry* hole = GEOSGetInteriorRingN(f->geom, h);
1171 "Looking for hole %d/%d of face %d among %d other faces",
1176 for (j = i + 1; j < nfaces; ++j)
1178 const GEOSGeometry* f2er;
1179 Face* f2 = faces[j];
1180 if (f2->parent)
continue;
1181 f2er = GEOSGetExteriorRing(f2->geom);
1184 if (GEOSEquals(f2er, hole))
1186 LWDEBUGF(2,
"Hole %d/%d of face %d is face %d", h + 1, nholes, i, j);
1195 static GEOSGeometry*
1196 collectFacesWithEvenAncestors(Face** faces,
int nfaces)
1198 GEOSGeometry** geoms =
lwalloc(
sizeof(GEOSGeometry*) * nfaces);
1200 unsigned int ngeoms = 0;
1203 for (i = 0; i < nfaces; ++i)
1206 if (countParens(f) % 2)
continue;
1207 geoms[ngeoms++] = GEOSGeom_clone(f->geom);
1210 ret = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, ngeoms);
1216 LWGEOM_GEOS_buildArea(
const GEOSGeometry* geom_in)
1219 GEOSGeometry *geos_result, *shp;
1220 GEOSGeometry
const* vgeoms[1];
1222 int srid = GEOSGetSRID(geom_in);
1224 #if POSTGIS_DEBUG_LEVEL >= 3
1229 vgeoms[0] = geom_in;
1230 geos_result = GEOSPolygonize(vgeoms, 1);
1232 LWDEBUGF(3,
"GEOSpolygonize returned @ %p", geos_result);
1235 if (!geos_result)
return 0;
1238 #if PARANOIA_LEVEL > 0
1241 GEOSGeom_destroy(geos_result);
1242 lwerror(
"%s [%d] Unexpected return from GEOSpolygonize", __FILE__, __LINE__);
1247 ngeoms = GEOSGetNumGeometries(geos_result);
1249 #if POSTGIS_DEBUG_LEVEL >= 3
1250 LWDEBUGF(3,
"GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
1253 LWDEBUGF(3,
"GEOSpolygonize: polygonized:%s", geom_ewkt);
1261 GEOSSetSRID(geos_result, srid);
1268 tmp = (GEOSGeometry*)GEOSGetGeometryN(geos_result, 0);
1271 GEOSGeom_destroy(geos_result);
1274 shp = GEOSGeom_clone(tmp);
1275 GEOSGeom_destroy(geos_result);
1276 GEOSSetSRID(shp, srid);
1280 LWDEBUGF(2,
"Polygonize returned %d geoms", ngeoms);
1310 geoms =
lwalloc(
sizeof(Face**) * ngeoms);
1311 for (i = 0; i < ngeoms; ++i)
1312 geoms[i] = newFace(GEOSGetGeometryN(geos_result, i));
1315 findFaceHoles(geoms, ngeoms);
1318 tmp = collectFacesWithEvenAncestors(geoms, ngeoms);
1321 for (i = 0; i < ngeoms; ++i)
1326 GEOSGeom_destroy(geos_result);
1329 shp = GEOSUnionCascaded(tmp);
1332 GEOSGeom_destroy(tmp);
1336 GEOSGeom_destroy(tmp);
1338 GEOSSetSRID(shp, srid);
1350 GEOSGeometry *g1, *g3;
1361 #if POSTGIS_GEOS_VERSION < 38
1362 g3 = LWGEOM_GEOS_buildArea(g1);
1364 g3 = GEOSBuildArea(g1);
1368 GEOSSetSRID(g3, srid);
1371 if (GEOSGetNumGeometries(g3) == 0)
1400 simple = GEOSisSimple(g);
1401 GEOSGeom_destroy(g);
1427 GEOSSetSRID(g, srid);
1443 GEOSGeometry *g1, *g2, *g3;
1452 g3 = GEOSSnap(g1, g2, tolerance);
1455 GEOSSetSRID(g3, srid);
1470 GEOSGeometry *g1, *g2, *g3;
1479 g3 = GEOSSharedPaths(g1, g2);
1482 GEOSSetSRID(g3, srid);
1498 GEOSGeometry *g1, *g3;
1506 g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
1514 GEOSSetSRID(g3, srid);
1536 for (i = 0; i < col->
ngeoms; i++)
1611 lwerror(
"lwgeom_offsetcurve: cannot node input");
1619 lwerror(
"lwgeom_offsetcurve: noded geometry cannot be offset");
1630 double area, bbox_area, bbox_width, bbox_height;
1633 uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
1634 double sample_cell_size;
1636 uint32_t iterations = 0;
1637 uint32_t npoints_generated = 0;
1638 uint32_t npoints_tested = 0;
1640 const GEOSPreparedGeometry* gprep;
1642 GEOSCoordSequence* gseq;
1647 const size_t size = 2 *
sizeof(int);
1648 char tmp[2 *
sizeof(int)];
1649 const size_t stride = 2 *
sizeof(int);
1653 lwerror(
"%s: only polygons supported", __func__);
1662 bbox = *(lwpoly->
bbox);
1665 bbox_width = bbox.
xmax - bbox.
xmin;
1666 bbox_height = bbox.
ymax - bbox.
ymin;
1667 bbox_area = bbox_width * bbox_height;
1669 if (area == 0.0 || bbox_area == 0.0)
1671 lwerror(
"%s: zero area input polygon, TBD", __func__);
1676 sample_npoints = npoints * bbox_area / area;
1681 sample_sqrt = lround(sqrt(sample_npoints));
1682 if (sample_sqrt == 0)
1686 if (bbox_width > bbox_height)
1688 sample_width = sample_sqrt;
1689 sample_height = ceil((
double)sample_npoints / (
double)sample_width);
1690 sample_cell_size = bbox_width / sample_width;
1694 sample_height = sample_sqrt;
1695 sample_width = ceil((
double)sample_npoints / (
double)sample_height);
1696 sample_cell_size = bbox_height / sample_height;
1707 gprep = GEOSPrepare(g);
1721 cells =
lwalloc(2 *
sizeof(
int) * sample_height * sample_width);
1722 for (i = 0; i < sample_width; i++)
1724 for (j = 0; j < sample_height; j++)
1726 cells[2 * (i * sample_height + j)] = i;
1727 cells[2 * (i * sample_height + j) + 1] = j;
1732 n = sample_height * sample_width;
1735 for (i = n - 1; i > 0; i--)
1739 memcpy(tmp, (
char *)cells + j * stride, size);
1740 memcpy((
char *)cells + j * stride, (
char *)cells + i * stride, size);
1741 memcpy((
char *)cells + i * stride, tmp, size);
1746 while (npoints_generated < npoints)
1749 for (i = 0; i < sample_width * sample_height; i++)
1752 double y = bbox.
ymin + cells[2 * i] * sample_cell_size;
1753 double x = bbox.
xmin + cells[2 * i + 1] * sample_cell_size;
1756 if (
x >= bbox.
xmax ||
y >= bbox.
ymax)
continue;
1758 gseq = GEOSCoordSeq_create(1, 2);
1759 #if POSTGIS_GEOS_VERSION < 38
1760 GEOSCoordSeq_setX(gseq, 0,
x);
1761 GEOSCoordSeq_setY(gseq, 0,
y);
1763 GEOSCoordSeq_setXY(gseq, 0,
x,
y);
1765 gpt = GEOSGeom_createPoint(gseq);
1767 contains = GEOSPreparedIntersects(gprep, gpt);
1769 GEOSGeom_destroy(gpt);
1773 GEOSPreparedGeom_destroy(gprep);
1774 GEOSGeom_destroy(g);
1780 npoints_generated++;
1782 if (npoints_generated == npoints)
1791 if (npoints_tested % 10000 == 0)
1792 LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g);
return NULL);
1796 if (done || iterations > 100)
break;
1799 GEOSPreparedGeom_destroy(gprep);
1800 GEOSGeom_destroy(g);
1818 lwerror(
"%s: only multipolygons supported", __func__);
1825 for (i = 0; i < lwmpoly->
ngeoms; i++)
1828 int sub_npoints = lround(npoints * sub_area / area);
1829 if (sub_npoints > 0)
1837 for (j = 0; j < sub_mpt->
ngeoms; j++)
1867 int type = GEOSGeomTypeId(geom);
1868 int SRID = GEOSGetSRID(geom);
1873 if (want3d && !GEOSHasZ(geom))
1875 LWDEBUG(3,
"Geometry has no Z, won't provide one");
1883 case GEOS_GEOMETRYCOLLECTION:
1884 LWDEBUG(4,
"lwgeom_from_geometry: it's a Collection or Multi");
1886 ngeoms = GEOSGetNumGeometries(geom);
1890 geoms =
lwalloc(ngeoms *
sizeof *geoms);
1893 lwerror(
"lwtin_from_geos: can't allocate geoms");
1896 for (i = 0; i < ngeoms; i++)
1898 const GEOSGeometry *poly, *ring;
1899 const GEOSCoordSequence* cs;
1902 poly = GEOSGetGeometryN(geom, i);
1903 ring = GEOSGetExteriorRing(poly);
1904 cs = GEOSGeom_getCoordSeq(ring);
1912 case GEOS_MULTIPOINT:
1913 case GEOS_MULTILINESTRING:
1914 case GEOS_MULTIPOLYGON:
1915 case GEOS_LINESTRING:
1916 case GEOS_LINEARRING:
1918 lwerror(
"lwtin_from_geos: invalid geometry type for tin: %d",
type);
1922 lwerror(
"GEOS2LWGEOM: unknown geometry type: %d",
type);
1938 GEOSGeometry *g1, *g3;
1940 if (output < 0 || output > 2)
1942 lwerror(
"%s: invalid output type specified %d", __func__, output);
1953 g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
1956 GEOSSetSRID(g3, srid);
1964 lwerror(
"%s: cannot convert output geometry", __func__);
1976 static GEOSCoordSequence*
1980 uint8_t num_dims = 2;
1982 GEOSCoordSequence* coords;
1985 coords = GEOSCoordSeq_create(num_points, num_dims);
1986 if (!coords)
return NULL;
1991 if (i >= num_points)
1993 lwerror(
"Incorrect num_points provided to lwgeom_get_geos_coordseq_2d");
1994 GEOSCoordSeq_destroy(coords);
1999 #if POSTGIS_GEOS_VERSION < 38
2000 if (!GEOSCoordSeq_setX(coords, i, tmp.
x) || !GEOSCoordSeq_setY(coords, i, tmp.
y))
2002 if (!GEOSCoordSeq_setXY(coords, i, tmp.
x, tmp.
y))
2005 GEOSCoordSeq_destroy(coords);
2023 GEOSCoordSequence* coords;
2024 GEOSGeometry* geos_geom;
2025 GEOSGeometry* geos_env = NULL;
2026 GEOSGeometry* geos_result;
2041 if (!coords)
return NULL;
2043 geos_geom = GEOSGeom_createLineString(coords);
2046 GEOSCoordSeq_destroy(coords);
2052 geos_result = GEOSVoronoiDiagram(geos_geom, geos_env, tolerance, output_edges);
2054 GEOSGeom_destroy(geos_geom);
2055 if (env) GEOSGeom_destroy(geos_env);
2064 GEOSGeom_destroy(geos_result);
2068 return lwgeom_result;
char result[OUT_DOUBLE_BUFFER_SIZE]
#define LWGEOM_GEOS_ERRMSG_MAXSIZE
LWGEOM * lwgeom_centroid(const LWGEOM *geom)
LWGEOM * lwgeom_intersection_prec(const LWGEOM *geom1, const LWGEOM *geom2, double prec)
LWGEOM * lwgeom_geos_noop(const LWGEOM *geom)
Convert an LWGEOM to a GEOS Geometry and convert back – for debug only.
static LWGEOM * lwcollection_offsetcurve(const LWCOLLECTION *col, double size, int quadsegs, int joinStyle, double mitreLimit)
static GEOSGeometry * ptarray_to_GEOSLinearRing(const POINTARRAY *pa, uint8_t autofix)
LWGEOM * lwgeom_symdifference(const LWGEOM *geom1, const LWGEOM *geom2)
LWGEOM * lwgeom_unaryunion(const LWGEOM *geom)
LWGEOM * lwgeom_union(const LWGEOM *g1, const LWGEOM *g2)
GEOSGeometry * make_geos_point(double x, double y)
static LWGEOM * lwline_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
LWGEOM * lwgeom_offsetcurve(const LWGEOM *geom, double size, int quadsegs, int joinStyle, double mitreLimit)
LWGEOM * lwgeom_delaunay_triangulation(const LWGEOM *geom, double tolerance, int32_t output)
Take vertices of a geometry and build a delaunay triangulation on them.
LWGEOM * lwgeom_difference_prec(const LWGEOM *geom1, const LWGEOM *geom2, double prec)
LWGEOM * lwgeom_snap(const LWGEOM *geom1, const LWGEOM *geom2, double tolerance)
Snap vertices and segments of a geometry to another using a given tolerance.
LWGEOM * lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2)
int lwgeom_is_simple(const LWGEOM *geom)
LWGEOM * lwgeom_pointonsurface(const LWGEOM *geom)
LWGEOM * lwgeom_sharedpaths(const LWGEOM *geom1, const LWGEOM *geom2)
LWGEOM * lwgeom_symdifference_prec(const LWGEOM *geom1, const LWGEOM *geom2, double prec)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
LWGEOM * lwgeom_voronoi_diagram(const LWGEOM *g, const GBOX *env, double tolerance, int output_edges)
Take vertices of a geometry and build the Voronoi diagram.
LWGEOM * lwgeom_union_prec(const LWGEOM *geom1, const LWGEOM *geom2, double gridSize)
LWGEOM * lwgeom_normalize(const LWGEOM *geom)
static int32_t get_result_srid(size_t count, const char *funcname,...)
#define GEOS_FREE_AND_FAIL_DEBUG(...)
LWGEOM * lwgeom_buildarea(const LWGEOM *geom)
Take a geometry and return an areal geometry (Polygon or MultiPolygon).
#define GEOS_FAIL_DEBUG()
LWGEOM * lwgeom_linemerge(const LWGEOM *geom)
#define GEOS_FREE_AND_FAIL(...)
GEOSGeometry * make_geos_segment(double x1, double y1, double x2, double y2)
static void geos_destroy(size_t count,...)
GEOSGeometry * GBOX2GEOS(const GBOX *box)
const char * lwgeom_geos_version()
Return GEOS version string (not to be freed)
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, uint8_t autofix)
LWGEOM * lwgeom_reduceprecision(const LWGEOM *geom, double gridSize)
LWMPOINT * lwgeom_to_points(const LWGEOM *lwgeom, uint32_t npoints, int32_t seed)
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, uint8_t want3d)
GEOSCoordSeq ptarray_to_GEOSCoordSeq(const POINTARRAY *, uint8_t fix_ring)
LWMPOINT * lwmpoly_to_points(const LWMPOLY *lwmpoly, uint32_t npoints, int32_t seed)
LWGEOM * lwgeom_unaryunion_prec(const LWGEOM *geom, double prec)
POINTARRAY * ptarray_from_GEOSCoordSeq(const GEOSCoordSequence *cs, uint8_t want3d)
LWMPOINT * lwpoly_to_points(const LWPOLY *lwpoly, uint32_t npoints, int32_t seed)
LWGEOM * lwgeom_intersection(const LWGEOM *g1, const LWGEOM *g2)
LWGEOM * lwgeom_clip_by_rect(const LWGEOM *geom1, double x1, double y1, double x2, double y2)
void lwgeom_geos_error(const char *fmt,...)
static GEOSCoordSequence * lwgeom_get_geos_coordseq_2d(const LWGEOM *g, uint32_t num_points)
LWTIN * lwtin_from_geos(const GEOSGeometry *geom, uint8_t want3d)
LWPOINT * lwpoint_construct_empty(int32_t srid, char hasz, char hasm)
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
LWGEOM * lwgeom_node(const LWGEOM *lwgeom_in)
LWPOINT * lwpoint_make2d(int32_t srid, double x, double y)
LWPOINTITERATOR * lwpointiterator_create(const LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM*.
void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
void lwgeom_free(LWGEOM *geom)
int lwpointiterator_next(LWPOINTITERATOR *s, POINT4D *p)
Attempts to assign the next point in the iterator to p, and advances the iterator to the next point.
POINTARRAY * ptarray_addPoint(const POINTARRAY *pa, uint8_t *p, size_t pdims, uint32_t where)
Add a point in a pointarray.
LWGEOM * lwgeom_stroke(const LWGEOM *geom, uint32_t perQuad)
Convert type with arcs into equivalent linearized type.
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
void lwpointiterator_destroy(LWPOINTITERATOR *s)
Free all memory associated with the iterator.
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
#define FLAGS_GET_Z(flags)
LWLINE * lwline_construct(int32_t srid, GBOX *bbox, POINTARRAY *points)
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
double lwgeom_area(const LWGEOM *geom)
int lwgeom_type_arc(const LWGEOM *geom)
Geometry type is one of the potentially "arc containing" types (circstring, multicurve,...
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
int lwgeom_dimension(const LWGEOM *geom)
For an LWGEOM, returns 0 for points, 1 for lines, 2 for polygons, 3 for volume, and the max dimension...
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
#define FLAGS_NDIMS(flags)
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
LWMPOINT * lwmpoint_construct_empty(int32_t srid, char hasz, char hasm)
LWPOINT * lwpoint_construct(int32_t srid, GBOX *bbox, POINTARRAY *point)
void lwcollection_release(LWCOLLECTION *lwcollection)
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int32_t srid, char hasz, char hasm)
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
void lwcollection_free(LWCOLLECTION *col)
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
LWCOLLECTION * lwcollection_concat_in_place(LWCOLLECTION *col1, const LWCOLLECTION *col2)
Appends all geometries from col2 to col1 in place.
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
void * lwalloc(size_t size)
LWCOLLECTION * lwcollection_construct(uint8_t type, int32_t srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
int ptarray_is_closed_2d(const POINTARRAY *pa)
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
#define LW_TRUE
Return types for functions with status returns.
void lwgeom_release(LWGEOM *lwgeom)
Free the containing LWGEOM and the associated BOX.
#define SRID_UNKNOWN
Unknown SRID value.
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
LWPOLY * lwpoly_construct_empty(int32_t srid, char hasz, char hasm)
LWPOLY * lwpoly_construct(int32_t srid, GBOX *bbox, uint32_t nrings, POINTARRAY **points)
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
LWLINE * lwline_construct_empty(int32_t srid, char hasz, char hasm)
LWTRIANGLE * lwtriangle_construct(int32_t srid, GBOX *bbox, POINTARRAY *points)
This library is the generic geometry handling section of PostGIS.
#define LW_ON_INTERRUPT(x)
double lwpoly_area(const LWPOLY *poly)
Find the area of the outer ring - sum (area of inner rings).
#define LWDEBUG(level, msg)
#define LWDEBUGF(level, msg,...)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
static const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from.
static const POINT3D * getPoint3d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT3D pointer into the POINTARRAY serialized_ptlist, suitable for reading from.
static uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
static uint8_t * getPoint_internal(const POINTARRAY *pa, uint32_t n)
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
void lwrandom_set_seed(int32_t seed)
double lwrandom_uniform(void)
Datum contains(PG_FUNCTION_ARGS)