PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwgeom_geos_split.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: lwgeom_geos.c 5258 2010-02-17 21:02:49Z strk $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  *
7  * Copyright 2009-2010 Sandro Santilli <strk@keybit.net>
8  *
9  * This is free software; you can redistribute and/or modify it under
10  * the terms of the GNU General Public Licence. See the COPYING file.
11  *
12  **********************************************************************
13  *
14  * Split polygon by line, line by line, line by point.
15  * Returns at most components as a collection.
16  * First element of the collection is always the part which
17  * remains after the cut, while the second element is the
18  * part which has been cut out. We arbitrarely take the part
19  * on the *right* of cut lines as the part which has been cut out.
20  * For a line cut by a point the part which remains is the one
21  * from start of the line to the cut point.
22  *
23  *
24  * Author: Sandro Santilli <strk@keybit.net>
25  *
26  * Work done for Faunalia (http://www.faunalia.it) with fundings
27  * from Regione Toscana - Sistema Informativo per il Governo
28  * del Territorio e dell'Ambiente (RT-SIGTA).
29  *
30  * Thanks to the PostGIS community for sharing poly/line ideas [1]
31  *
32  * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
33  *
34  *
35  **********************************************************************/
36 
37 #include "lwgeom_geos.h"
38 #include "liblwgeom_internal.h"
39 
40 #include <string.h>
41 #include <assert.h>
42 
43 static LWGEOM* lwline_split_by_line(const LWLINE* lwgeom_in, const LWLINE* blade_in);
44 static LWGEOM* lwline_split_by_point(const LWLINE* lwgeom_in, const LWPOINT* blade_in);
45 static LWGEOM* lwline_split(const LWLINE* lwgeom_in, const LWGEOM* blade_in);
46 static LWGEOM* lwpoly_split_by_line(const LWPOLY* lwgeom_in, const LWLINE* blade_in);
47 static LWGEOM* lwcollection_split(const LWCOLLECTION* lwcoll_in, const LWGEOM* blade_in);
48 static LWGEOM* lwpoly_split(const LWPOLY* lwpoly_in, const LWGEOM* blade_in);
49 
50 /* Initializes and uses GEOS internally */
51 static LWGEOM*
52 lwline_split_by_line(const LWLINE* lwline_in, const LWLINE* blade_in)
53 {
54  LWGEOM** components;
55  LWGEOM* diff;
56  LWCOLLECTION* out;
57  GEOSGeometry* gdiff; /* difference */
58  GEOSGeometry* g1;
59  GEOSGeometry* g2;
60  int ret;
61 
62  /* Possible outcomes:
63  *
64  * 1. The lines do not cross or overlap
65  * -> Return a collection with single element
66  * 2. The lines cross
67  * -> Return a collection of all elements resulting from the split
68  */
69 
71 
72  g1 = LWGEOM2GEOS((LWGEOM*)lwline_in);
73  if ( ! g1 )
74  {
75  lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
76  return NULL;
77  }
78  g2 = LWGEOM2GEOS((LWGEOM*)blade_in);
79  if ( ! g2 )
80  {
81  GEOSGeom_destroy(g1);
82  lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
83  return NULL;
84  }
85 
86  /* If interior intersecton is linear we can't split */
87  ret = GEOSRelatePattern(g1, g2, "1********");
88  if ( 2 == ret )
89  {
90  lwerror("GEOSRelatePattern: %s", lwgeom_geos_errmsg);
91  GEOSGeom_destroy(g1);
92  GEOSGeom_destroy(g2);
93  return NULL;
94  }
95  if ( ret )
96  {
97  GEOSGeom_destroy(g1);
98  GEOSGeom_destroy(g2);
99  lwerror("Splitter line has linear intersection with input");
100  return NULL;
101  }
102 
103 
104  gdiff = GEOSDifference(g1,g2);
105  GEOSGeom_destroy(g1);
106  GEOSGeom_destroy(g2);
107  if (gdiff == NULL)
108  {
109  lwerror("GEOSDifference: %s", lwgeom_geos_errmsg);
110  return NULL;
111  }
112 
113  diff = GEOS2LWGEOM(gdiff, FLAGS_GET_Z(lwline_in->flags));
114  GEOSGeom_destroy(gdiff);
115  if (NULL == diff)
116  {
117  lwerror("GEOS2LWGEOM: %s", lwgeom_geos_errmsg);
118  return NULL;
119  }
120 
121  out = lwgeom_as_lwcollection(diff);
122  if ( ! out )
123  {
124  components = lwalloc(sizeof(LWGEOM*)*1);
125  components[0] = diff;
126  out = lwcollection_construct(COLLECTIONTYPE, lwline_in->srid,
127  NULL, 1, components);
128  }
129  else
130  {
131  /* Set SRID */
132  lwgeom_set_srid((LWGEOM*)out, lwline_in->srid);
133  /* Force collection type */
134  out->type = COLLECTIONTYPE;
135  }
136 
137 
138  return (LWGEOM*)out;
139 }
140 
141 static LWGEOM*
142 lwline_split_by_point(const LWLINE* lwline_in, const LWPOINT* blade_in)
143 {
144  LWMLINE* out;
145 
146  out = lwmline_construct_empty(lwline_in->srid,
147  FLAGS_GET_Z(lwline_in->flags),
148  FLAGS_GET_M(lwline_in->flags));
149  if ( lwline_split_by_point_to(lwline_in, blade_in, out) < 2 )
150  {
151  lwmline_add_lwline(out, lwline_clone_deep(lwline_in));
152  }
153 
154  /* Turn multiline into collection */
155  out->type = COLLECTIONTYPE;
156 
157  return (LWGEOM*)out;
158 }
159 
160 int
161 lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
162  LWMLINE* v)
163 {
164  double loc, dist;
165  POINT4D pt, pt_projected;
166  POINTARRAY* pa1;
167  POINTARRAY* pa2;
168  double vstol; /* vertex snap tolerance */
169 
170  /* Possible outcomes:
171  *
172  * 1. The point is not on the line or on the boundary
173  * -> Leave collection untouched, return 0
174  * 2. The point is on the boundary
175  * -> Leave collection untouched, return 1
176  * 3. The point is in the line
177  * -> Push 2 elements on the collection:
178  * o start_point - cut_point
179  * o cut_point - last_point
180  * -> Return 2
181  */
182 
183  getPoint4d_p(blade_in->point, 0, &pt);
184  loc = ptarray_locate_point(lwline_in->points, &pt, &dist, &pt_projected);
185 
186  /* lwnotice("Location: %g -- Distance: %g", loc, dist); */
187 
188  if ( dist > 0 ) /* TODO: accept a tolerance ? */
189  {
190  /* No intersection */
191  return 0;
192  }
193 
194  if ( loc == 0 || loc == 1 )
195  {
196  /* Intersection is on the boundary */
197  return 1;
198  }
199 
200  /* There is a real intersection, let's get two substrings */
201 
202  /* Compute vertex snap tolerance based on line length
203  * TODO: take as parameter ? */
204  vstol = ptarray_length_2d(lwline_in->points) / 1e14;
205 
206  pa1 = ptarray_substring(lwline_in->points, 0, loc, vstol);
207  pa2 = ptarray_substring(lwline_in->points, loc, 1, vstol);
208 
209  /* NOTE: I've seen empty pointarrays with loc != 0 and loc != 1 */
210  if ( pa1->npoints == 0 || pa2->npoints == 0 ) {
211  ptarray_free(pa1);
212  ptarray_free(pa2);
213  /* Intersection is on the boundary */
214  return 1;
215  }
216 
219  return 2;
220 }
221 
222 static LWGEOM*
223 lwline_split(const LWLINE* lwline_in, const LWGEOM* blade_in)
224 {
225  switch (blade_in->type)
226  {
227  case POINTTYPE:
228  return lwline_split_by_point(lwline_in, (LWPOINT*)blade_in);
229 
230  case LINETYPE:
231  return lwline_split_by_line(lwline_in, (LWLINE*)blade_in);
232 
233  default:
234  lwerror("Splitting a Line by a %s is unsupported",
235  lwtype_name(blade_in->type));
236  return NULL;
237  }
238  return NULL;
239 }
240 
241 /* Initializes and uses GEOS internally */
242 static LWGEOM*
243 lwpoly_split_by_line(const LWPOLY* lwpoly_in, const LWLINE* blade_in)
244 {
245  LWCOLLECTION* out;
246  GEOSGeometry* g1;
247  GEOSGeometry* g2;
248  GEOSGeometry* g1_bounds;
249  GEOSGeometry* polygons;
250  const GEOSGeometry *vgeoms[1];
251  int i,n;
252  int hasZ = FLAGS_GET_Z(lwpoly_in->flags);
253 
254 
255  /* Possible outcomes:
256  *
257  * 1. The line does not split the polygon
258  * -> Return a collection with single element
259  * 2. The line does split the polygon
260  * -> Return a collection of all elements resulting from the split
261  */
262 
264 
265  g1 = LWGEOM2GEOS((LWGEOM*)lwpoly_in);
266  if ( NULL == g1 )
267  {
268  lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
269  return NULL;
270  }
271  g1_bounds = GEOSBoundary(g1);
272  if ( NULL == g1_bounds )
273  {
274  GEOSGeom_destroy(g1);
275  lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg);
276  return NULL;
277  }
278 
279  g2 = LWGEOM2GEOS((LWGEOM*)blade_in);
280  if ( NULL == g2 )
281  {
282  GEOSGeom_destroy(g1);
283  GEOSGeom_destroy(g1_bounds);
284  lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
285  return NULL;
286  }
287 
288  vgeoms[0] = GEOSUnion(g1_bounds, g2);
289  if ( NULL == vgeoms[0] )
290  {
291  GEOSGeom_destroy(g1);
292  GEOSGeom_destroy(g2);
293  GEOSGeom_destroy(g1_bounds);
294  lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
295  return NULL;
296  }
297 
298  /* debugging..
299  lwnotice("Bounds poly: %s",
300  lwgeom_to_ewkt(GEOS2LWGEOM(g1_bounds, hasZ)));
301  lwnotice("Line: %s",
302  lwgeom_to_ewkt(GEOS2LWGEOM(g2, hasZ)));
303 
304  lwnotice("Noded bounds: %s",
305  lwgeom_to_ewkt(GEOS2LWGEOM(vgeoms[0], hasZ)));
306  */
307 
308  polygons = GEOSPolygonize(vgeoms, 1);
309  if ( NULL == polygons )
310  {
311  GEOSGeom_destroy(g1);
312  GEOSGeom_destroy(g2);
313  GEOSGeom_destroy(g1_bounds);
314  GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
315  lwerror("GEOSPolygonize: %s", lwgeom_geos_errmsg);
316  return NULL;
317  }
318 
319 #if PARANOIA_LEVEL > 0
320  if ( GEOSGeometryTypeId(polygons) != COLLECTIONTYPE )
321  {
322  GEOSGeom_destroy(g1);
323  GEOSGeom_destroy(g2);
324  GEOSGeom_destroy(g1_bounds);
325  GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
326  GEOSGeom_destroy(polygons);
327  lwerror("Unexpected return from GEOSpolygonize");
328  return 0;
329  }
330 #endif
331 
332  /* We should now have all polygons, just skip
333  * the ones which are in holes of the original
334  * geometries and return the rest in a collection
335  */
336  n = GEOSGetNumGeometries(polygons);
338  hasZ, 0);
339  /* Allocate space for all polys */
340  out->geoms = lwrealloc(out->geoms, sizeof(LWGEOM*)*n);
341  assert(0 == out->ngeoms);
342  for (i=0; i<n; ++i)
343  {
344  GEOSGeometry* pos; /* point on surface */
345  const GEOSGeometry* p = GEOSGetGeometryN(polygons, i);
346  int contains;
347 
348  pos = GEOSPointOnSurface(p);
349  if ( ! pos )
350  {
351  GEOSGeom_destroy(g1);
352  GEOSGeom_destroy(g2);
353  GEOSGeom_destroy(g1_bounds);
354  GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
355  GEOSGeom_destroy(polygons);
356  lwerror("GEOSPointOnSurface: %s", lwgeom_geos_errmsg);
357  return NULL;
358  }
359 
360  contains = GEOSContains(g1, pos);
361  if ( 2 == contains )
362  {
363  GEOSGeom_destroy(g1);
364  GEOSGeom_destroy(g2);
365  GEOSGeom_destroy(g1_bounds);
366  GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
367  GEOSGeom_destroy(polygons);
368  GEOSGeom_destroy(pos);
369  lwerror("GEOSContains: %s", lwgeom_geos_errmsg);
370  return NULL;
371  }
372 
373  GEOSGeom_destroy(pos);
374 
375  if ( 0 == contains )
376  {
377  /* Original geometry doesn't contain
378  * a point in this ring, must be an hole
379  */
380  continue;
381  }
382 
383  out->geoms[out->ngeoms++] = GEOS2LWGEOM(p, hasZ);
384  }
385 
386  GEOSGeom_destroy(g1);
387  GEOSGeom_destroy(g2);
388  GEOSGeom_destroy(g1_bounds);
389  GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
390  GEOSGeom_destroy(polygons);
391 
392  return (LWGEOM*)out;
393 }
394 
395 static LWGEOM*
396 lwcollection_split(const LWCOLLECTION* lwcoll_in, const LWGEOM* blade_in)
397 {
398  LWGEOM** split_vector=NULL;
399  LWCOLLECTION* out;
400  size_t split_vector_capacity;
401  size_t split_vector_size=0;
402  size_t i,j;
403 
404  split_vector_capacity=8;
405  split_vector = lwalloc(split_vector_capacity * sizeof(LWGEOM*));
406  if ( ! split_vector )
407  {
408  lwerror("Out of virtual memory");
409  return NULL;
410  }
411 
412  for (i=0; i<lwcoll_in->ngeoms; ++i)
413  {
414  LWCOLLECTION* col;
415  LWGEOM* split = lwgeom_split(lwcoll_in->geoms[i], blade_in);
416  /* an exception should prevent this from ever returning NULL */
417  if ( ! split ) return NULL;
418 
419  col = lwgeom_as_lwcollection(split);
420  /* Output, if any, will always be a collection */
421  assert(col);
422 
423  /* Reallocate split_vector if needed */
424  if ( split_vector_size + col->ngeoms > split_vector_capacity )
425  {
426  /* NOTE: we could be smarter on reallocations here */
427  split_vector_capacity += col->ngeoms;
428  split_vector = lwrealloc(split_vector,
429  split_vector_capacity * sizeof(LWGEOM*));
430  if ( ! split_vector )
431  {
432  lwerror("Out of virtual memory");
433  return NULL;
434  }
435  }
436 
437  for (j=0; j<col->ngeoms; ++j)
438  {
439  col->geoms[j]->srid = SRID_UNKNOWN; /* strip srid */
440  split_vector[split_vector_size++] = col->geoms[j];
441  }
442  lwfree(col->geoms);
443  lwfree(col);
444  }
445 
446  /* Now split_vector has split_vector_size geometries */
447  out = lwcollection_construct(COLLECTIONTYPE, lwcoll_in->srid,
448  NULL, split_vector_size, split_vector);
449 
450  return (LWGEOM*)out;
451 }
452 
453 static LWGEOM*
454 lwpoly_split(const LWPOLY* lwpoly_in, const LWGEOM* blade_in)
455 {
456  switch (blade_in->type)
457  {
458  case LINETYPE:
459  return lwpoly_split_by_line(lwpoly_in, (LWLINE*)blade_in);
460  default:
461  lwerror("Splitting a Polygon by a %s is unsupported",
462  lwtype_name(blade_in->type));
463  return NULL;
464  }
465  return NULL;
466 }
467 
468 /* exported */
469 LWGEOM*
470 lwgeom_split(const LWGEOM* lwgeom_in, const LWGEOM* blade_in)
471 {
472  switch (lwgeom_in->type)
473  {
474  case LINETYPE:
475  return lwline_split((const LWLINE*)lwgeom_in, blade_in);
476 
477  case POLYGONTYPE:
478  return lwpoly_split((const LWPOLY*)lwgeom_in, blade_in);
479 
480  case MULTIPOLYGONTYPE:
481  case MULTILINETYPE:
482  case COLLECTIONTYPE:
483  return lwcollection_split((const LWCOLLECTION*)lwgeom_in, blade_in);
484 
485  default:
486  lwerror("Splitting of %s geometries is unsupported",
487  lwtype_name(lwgeom_in->type));
488  return NULL;
489  }
490 
491 }
492 
#define LINETYPE
Definition: liblwgeom.h:61
LWGEOM * lwgeom_split(const LWGEOM *lwgeom_in, const LWGEOM *blade_in)
uint8_t type
Definition: liblwgeom.h:433
static LWGEOM * lwline_split(const LWLINE *lwgeom_in, const LWGEOM *blade_in)
LWCOLLECTION * lwcollection_construct(uint8_t type, int srid, GBOX *bbox, uint32_t ngeoms, LWGEOM **geoms)
Definition: lwcollection.c:30
static LWGEOM * lwline_split_by_line(const LWLINE *lwgeom_in, const LWLINE *blade_in)
void lwfree(void *mem)
Definition: lwutil.c:190
int npoints
Definition: liblwgeom.h:327
uint8_t type
Definition: liblwgeom.h:459
#define POLYGONTYPE
Definition: liblwgeom.h:62
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:315
POINTARRAY * ptarray_substring(POINTARRAY *pa, double d1, double d2, double tolerance)
start location (distance from start / total distance) end location (distance from start / total dist...
Definition: ptarray.c:1022
double ptarray_length_2d(const POINTARRAY *pts)
Find the 2d length of the given POINTARRAY (even if it's 3d)
Definition: ptarray.c:1586
Datum contains(PG_FUNCTION_ARGS)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom)
LWMLINE * lwmline_add_lwline(LWMLINE *mobj, const LWLINE *obj)
Definition: lwmline.c:33
int32_t srid
Definition: liblwgeom.h:377
POINTARRAY * point
Definition: liblwgeom.h:367
LWMLINE * lwmline_construct_empty(int srid, char hasz, char hasm)
Definition: lwmline.c:25
int32_t srid
Definition: liblwgeom.h:355
static LWGEOM * lwcollection_split(const LWCOLLECTION *lwcoll_in, const LWGEOM *blade_in)
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double ptarray_locate_point(const POINTARRAY *pa, const POINT4D *pt, double *dist, POINT4D *p_located)
Definition: ptarray.c:1267
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
void lwgeom_geos_error(const char *fmt,...)
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
LWGEOM ** geoms
Definition: liblwgeom.h:465
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:154
int32_t srid
Definition: liblwgeom.h:462
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
int lwline_split_by_point_to(const LWLINE *lwline_in, const LWPOINT *blade_in, LWMLINE *v)
Split a line by a point and push components to the provided multiline.
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
int32_t srid
Definition: liblwgeom.h:410
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:143
static LWGEOM * lwpoly_split(const LWPOLY *lwpoly_in, const LWGEOM *blade_in)
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:183
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:107
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
uint8_t type
Definition: liblwgeom.h:352
void lwgeom_set_srid(LWGEOM *geom, int srid)
Set the SRID on an LWGEOM For collections, only the parent gets an SRID, all the children get SRID_UN...
uint8_t flags
Definition: liblwgeom.h:408
void * lwalloc(size_t size)
Definition: lwutil.c:175
#define MULTILINETYPE
Definition: liblwgeom.h:64
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
LWLINE * lwline_clone_deep(const LWLINE *lwgeom)
Definition: lwline.c:105
uint8_t flags
Definition: liblwgeom.h:375
static LWGEOM * lwpoly_split_by_line(const LWPOLY *lwgeom_in, const LWLINE *blade_in)
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:217
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
static LWGEOM * lwline_split_by_point(const LWLINE *lwgeom_in, const LWPOINT *blade_in)
POINTARRAY * points
Definition: liblwgeom.h:378