PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwlinearreferencing.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id$
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  * Copyright 2011 Paul Ramsey
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "liblwgeom_internal.h"
14 #include "lwgeom_log.h"
15 
16 
17 static int
18 segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn)
19 {
20  double m1 = p1->m;
21  double m2 = p2->m;
22  double mprop;
23 
24  /* M is out of range, no new point generated. */
25  if ( (m < FP_MIN(m1,m2)) || (m > FP_MAX(m1,m2)) )
26  {
27  return LW_FALSE;
28  }
29 
30  /* We'll just can out on this degenerate case for now.
31  Correct behavior is probably an mprop of 0.5?
32  Still would have to deal with case of true p1==p2. */
33  if( m1 == m2 )
34  {
35  lwerror("Zero measure-length line encountered!");
36  }
37 
38  /* M is in range, new point to be generated. */
39  mprop = (m - m1) / (m2 - m1);
40  pn->x = p1->x + (p2->x - p1->x) * mprop;
41  pn->y = p1->y + (p2->y - p1->y) * mprop;
42  pn->z = p1->z + (p2->z - p1->z) * mprop;
43  pn->m = m;
44 
45  /* Offset to the left or right, if necessary. */
46  if ( offset != 0.0 )
47  {
48  double theta = atan2(p2->y - p1->y, p2->x - p1->x);
49  pn->x -= sin(theta) * offset;
50  pn->y += cos(theta) * offset;
51  }
52 
53  return LW_TRUE;
54 }
55 
56 
57 static POINTARRAY*
58 ptarray_locate_along(const POINTARRAY *pa, double m, double offset)
59 {
60  int i;
61  POINT4D p1, p2, pn;
62  POINTARRAY *dpa = NULL;
63 
64  /* Can't do anything with degenerate point arrays */
65  if ( ! pa || pa->npoints < 2 ) return NULL;
66 
67  /* Walk through each segment in the point array */
68  for ( i = 1; i < pa->npoints; i++ )
69  {
70  getPoint4d_p(pa, i-1, &p1);
71  getPoint4d_p(pa, i, &p2);
72 
73  /* No derived point? Move to next segment. */
74  if ( segment_locate_along(&p1, &p2, m, offset, &pn) == LW_FALSE )
75  continue;
76 
77  /* No pointarray, make a fresh one */
78  if ( dpa == NULL )
80 
81  /* Add our new point to the array */
82  ptarray_append_point(dpa, &pn, 0);
83  }
84 
85  return dpa;
86 }
87 
88 static LWMPOINT*
89 lwline_locate_along(const LWLINE *lwline, double m, double offset)
90 {
91  POINTARRAY *opa = NULL;
92  LWMPOINT *mp = NULL;
93  LWGEOM *lwg = lwline_as_lwgeom(lwline);
94  int hasz, hasm, srid;
95 
96  /* Return degenerates upwards */
97  if ( ! lwline ) return NULL;
98 
99  /* Create empty return shell */
100  srid = lwgeom_get_srid(lwg);
101  hasz = lwgeom_has_z(lwg);
102  hasm = lwgeom_has_m(lwg);
103 
104  if ( hasm )
105  {
106  /* Find points along */
107  opa = ptarray_locate_along(lwline->points, m, offset);
108  }
109  else
110  {
111  LWLINE *lwline_measured = lwline_measured_from_lwline(lwline, 0.0, 1.0);
112  opa = ptarray_locate_along(lwline_measured->points, m, offset);
113  lwline_free(lwline_measured);
114  }
115 
116  /* Return NULL as EMPTY */
117  if ( ! opa )
118  return lwmpoint_construct_empty(srid, hasz, hasm);
119 
120  /* Convert pointarray into a multipoint */
121  mp = lwmpoint_construct(srid, opa);
122  ptarray_free(opa);
123  return mp;
124 }
125 
126 static LWMPOINT*
127 lwmline_locate_along(const LWMLINE *lwmline, double m, double offset)
128 {
129  LWMPOINT *lwmpoint = NULL;
130  LWGEOM *lwg = lwmline_as_lwgeom(lwmline);
131  int i, j;
132 
133  /* Return degenerates upwards */
134  if ( (!lwmline) || (lwmline->ngeoms < 1) ) return NULL;
135 
136  /* Construct return */
138 
139  /* Locate along each sub-line */
140  for ( i = 0; i < lwmline->ngeoms; i++ )
141  {
142  LWMPOINT *along = lwline_locate_along(lwmline->geoms[i], m, offset);
143  if ( along )
144  {
145  if ( ! lwgeom_is_empty((LWGEOM*)along) )
146  {
147  for ( j = 0; j < along->ngeoms; j++ )
148  {
149  lwmpoint_add_lwpoint(lwmpoint, along->geoms[j]);
150  }
151  }
152  /* Free the containing geometry, but leave the sub-geometries around */
153  along->ngeoms = 0; lwmpoint_free(along);
154  }
155  }
156  return lwmpoint;
157 }
158 
159 static LWMPOINT*
160 lwpoint_locate_along(const LWPOINT *lwpoint, double m, double offset)
161 {
162  double point_m = lwpoint_get_m(lwpoint);
163  LWGEOM *lwg = lwpoint_as_lwgeom(lwpoint);
165  if ( FP_EQUALS(m, point_m) )
166  {
167  lwmpoint_add_lwpoint(r, lwpoint_clone(lwpoint));
168  }
169  return r;
170 }
171 
172 static LWMPOINT*
173 lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset)
174 {
175  LWGEOM *lwg = lwmpoint_as_lwgeom(lwin);
176  LWMPOINT *lwout = NULL;
177  int i;
178 
179  /* Construct return */
181 
182  for ( i = 0; i < lwin->ngeoms; i++ )
183  {
184  double point_m = lwpoint_get_m(lwin->geoms[i]);
185  if ( FP_EQUALS(m, point_m) )
186  {
187  lwmpoint_add_lwpoint(lwout, lwpoint_clone(lwin->geoms[i]));
188  }
189  }
190 
191  return lwout;
192 }
193 
194 LWGEOM*
195 lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
196 {
197  if ( ! lwin ) return NULL;
198 
199  if ( ! lwgeom_has_m(lwin) )
200  lwerror("Input geometry does not have a measure dimension");
201 
202  switch (lwin->type)
203  {
204  case POINTTYPE:
205  return (LWGEOM*)lwpoint_locate_along((LWPOINT*)lwin, m, offset);
206  case MULTIPOINTTYPE:
207  return (LWGEOM*)lwmpoint_locate_along((LWMPOINT*)lwin, m, offset);
208  case LINETYPE:
209  return (LWGEOM*)lwline_locate_along((LWLINE*)lwin, m, offset);
210  case MULTILINETYPE:
211  return (LWGEOM*)lwmline_locate_along((LWMLINE*)lwin, m, offset);
212  /* Only line types supported right now */
213  /* TO DO: CurveString, CompoundCurve, MultiCurve */
214  /* TO DO: Point, MultiPoint */
215  default:
216  lwerror("Only linear geometries are supported, %s provided.",lwtype_name(lwin->type));
217  return NULL;
218  }
219  return NULL;
220 }
221 
229 double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
230 {
231  if ( ! p )
232  {
233  lwerror("Null input geometry.");
234  return 0.0;
235  }
236 
237  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
238  {
239  lwerror("Cannot extract %c ordinate.", ordinate);
240  return 0.0;
241  }
242 
243  if ( ordinate == 'X' )
244  return p->x;
245  if ( ordinate == 'Y' )
246  return p->y;
247  if ( ordinate == 'Z' )
248  return p->z;
249  if ( ordinate == 'M' )
250  return p->m;
251 
252  /* X */
253  return p->x;
254 
255 }
256 
261 void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
262 {
263  if ( ! p )
264  {
265  lwerror("Null input geometry.");
266  return;
267  }
268 
269  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
270  {
271  lwerror("Cannot set %c ordinate.", ordinate);
272  return;
273  }
274 
275  LWDEBUGF(4, " setting ordinate %c to %g", ordinate, value);
276 
277  switch ( ordinate )
278  {
279  case 'X':
280  p->x = value;
281  return;
282  case 'Y':
283  p->y = value;
284  return;
285  case 'Z':
286  p->z = value;
287  return;
288  case 'M':
289  p->m = value;
290  return;
291  }
292 }
293 
299 int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
300 {
301  static char* dims = "XYZM";
302  double p1_value = lwpoint_get_ordinate(p1, ordinate);
303  double p2_value = lwpoint_get_ordinate(p2, ordinate);
304  double proportion;
305  int i = 0;
306 
307  if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
308  {
309  lwerror("Cannot set %c ordinate.", ordinate);
310  return 0;
311  }
312 
313  if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
314  FP_MAX(p1_value, p2_value) < interpolation_value )
315  {
316  lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
317  return 0;
318  }
319 
320  proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
321 
322  for ( i = 0; i < 4; i++ )
323  {
324  double newordinate = 0.0;
325  if ( dims[i] == 'Z' && ! hasz ) continue;
326  if ( dims[i] == 'M' && ! hasm ) continue;
327  p1_value = lwpoint_get_ordinate(p1, dims[i]);
328  p2_value = lwpoint_get_ordinate(p2, dims[i]);
329  newordinate = p1_value + proportion * (p2_value - p1_value);
330  lwpoint_set_ordinate(p, dims[i], newordinate);
331  LWDEBUGF(4, " clip ordinate(%c) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", dims[i], p1_value, p2_value, proportion, newordinate );
332  }
333 
334  return 1;
335 }
336 
337 
342 lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
343 {
344  LWCOLLECTION *lwgeom_out = NULL;
345  char hasz, hasm;
346  POINT4D p4d;
347  double ordinate_value;
348 
349  /* Nothing to do with NULL */
350  if ( ! point )
351  lwerror("Null input geometry.");
352 
353  /* Ensure 'from' is less than 'to'. */
354  if ( to < from )
355  {
356  double t = from;
357  from = to;
358  to = t;
359  }
360 
361  /* Read Z/M info */
362  hasz = lwgeom_has_z(lwpoint_as_lwgeom(point));
363  hasm = lwgeom_has_m(lwpoint_as_lwgeom(point));
364 
365  /* Prepare return object */
366  lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, point->srid, hasz, hasm);
367 
368  /* Test if ordinate is in range */
369  lwpoint_getPoint4d_p(point, &p4d);
370  ordinate_value = lwpoint_get_ordinate(&p4d, ordinate);
371  if ( from <= ordinate_value && to >= ordinate_value )
372  {
373  LWPOINT *lwp = lwpoint_clone(point);
374  lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
375  }
376 
377  /* Set the bbox, if necessary */
378  if ( lwgeom_out->bbox )
379  {
380  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
381  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
382  }
383 
384  return lwgeom_out;
385 }
386 
387 
388 
393 lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
394 {
395  LWCOLLECTION *lwgeom_out = NULL;
396  char hasz, hasm;
397  int i;
398 
399  /* Nothing to do with NULL */
400  if ( ! mpoint )
401  lwerror("Null input geometry.");
402 
403  /* Ensure 'from' is less than 'to'. */
404  if ( to < from )
405  {
406  double t = from;
407  from = to;
408  to = t;
409  }
410 
411  /* Read Z/M info */
412  hasz = lwgeom_has_z(lwmpoint_as_lwgeom(mpoint));
413  hasm = lwgeom_has_m(lwmpoint_as_lwgeom(mpoint));
414 
415  /* Prepare return object */
416  lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, mpoint->srid, hasz, hasm);
417 
418  /* For each point, is its ordinate value between from and to? */
419  for ( i = 0; i < mpoint->ngeoms; i ++ )
420  {
421  POINT4D p4d;
422  double ordinate_value;
423 
424  lwpoint_getPoint4d_p(mpoint->geoms[i], &p4d);
425  ordinate_value = lwpoint_get_ordinate(&p4d, ordinate);
426 
427  if ( from <= ordinate_value && to >= ordinate_value )
428  {
429  LWPOINT *lwp = lwpoint_clone(mpoint->geoms[i]);
430  lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
431  }
432  }
433 
434  /* Set the bbox, if necessary */
435  if ( lwgeom_out->bbox )
436  {
437  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
438  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
439  }
440 
441  return lwgeom_out;
442 }
443 
448 lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
449 {
450  LWCOLLECTION *lwgeom_out = NULL;
451 
452  if ( ! mline )
453  {
454  lwerror("Null input geometry.");
455  return NULL;
456  }
457 
458  if ( mline->ngeoms == 1)
459  {
460  lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
461  }
462  else
463  {
464  LWCOLLECTION *col;
465  char hasz = lwgeom_has_z(lwmline_as_lwgeom(mline));
466  char hasm = lwgeom_has_m(lwmline_as_lwgeom(mline));
467  int i, j;
468  char homogeneous = 1;
469  size_t geoms_size = 0;
470  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, mline->srid, hasz, hasm);
471  FLAGS_SET_Z(lwgeom_out->flags, hasz);
472  FLAGS_SET_M(lwgeom_out->flags, hasm);
473  for ( i = 0; i < mline->ngeoms; i ++ )
474  {
475  col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
476  if ( col )
477  {
478  /* Something was left after the clip. */
479  if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
480  {
481  geoms_size += 16;
482  if ( lwgeom_out->geoms )
483  {
484  lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
485  }
486  else
487  {
488  lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
489  }
490  }
491  for ( j = 0; j < col->ngeoms; j++ )
492  {
493  lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
494  lwgeom_out->ngeoms++;
495  }
496  if ( col->type != mline->type )
497  {
498  homogeneous = 0;
499  }
500  /* Shallow free the struct, leaving the geoms behind. */
501  if ( col->bbox ) lwfree(col->bbox);
502  lwfree(col->geoms);
503  lwfree(col);
504  }
505  }
506  if ( lwgeom_out->bbox )
507  {
508  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
509  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
510  }
511 
512  if ( ! homogeneous )
513  {
514  lwgeom_out->type = COLLECTIONTYPE;
515  }
516  }
517 
518  return lwgeom_out;
519 
520 }
521 
522 
528 lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
529 {
530 
531  POINTARRAY *pa_in = NULL;
532  LWCOLLECTION *lwgeom_out = NULL;
533  POINTARRAY *dp = NULL;
534  int i;
535  int added_last_point = 0;
536  POINT4D *p = NULL, *q = NULL, *r = NULL;
537  double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
538  char hasz = lwgeom_has_z(lwline_as_lwgeom(line));
539  char hasm = lwgeom_has_m(lwline_as_lwgeom(line));
540  char dims = FLAGS_NDIMS(line->flags);
541 
542  /* Null input, nothing we can do. */
543  if ( ! line )
544  {
545  lwerror("Null input geometry.");
546  return NULL;
547  }
548 
549  /* Ensure 'from' is less than 'to'. */
550  if ( to < from )
551  {
552  double t = from;
553  from = to;
554  to = t;
555  }
556 
557  LWDEBUGF(4, "from = %g, to = %g, ordinate = %c", from, to, ordinate);
558  LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
559 
560  /* Asking for an ordinate we don't have. Error. */
561  if ( (ordinate == 'Z' && ! hasz) || (ordinate == 'M' && ! hasm) )
562  {
563  lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
564  return NULL;
565  }
566 
567  /* Prepare our working point objects. */
568  p = lwalloc(sizeof(POINT4D));
569  q = lwalloc(sizeof(POINT4D));
570  r = lwalloc(sizeof(POINT4D));
571 
572  /* Construct a collection to hold our outputs. */
573  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
574 
575  /* Get our input point array */
576  pa_in = line->points;
577 
578  for ( i = 0; i < pa_in->npoints; i++ )
579  {
580  LWDEBUGF(4, "Point #%d", i);
581  LWDEBUGF(4, "added_last_point %d", added_last_point);
582  if ( i > 0 )
583  {
584  *q = *p;
585  ordinate_value_q = ordinate_value_p;
586  }
587  getPoint4d_p(pa_in, i, p);
588  ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
589  LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
590  LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
591 
592  /* Is this point inside the ordinate range? Yes. */
593  if ( ordinate_value_p >= from && ordinate_value_p <= to )
594  {
595  LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
596 
597  if ( ! added_last_point )
598  {
599  LWDEBUG(4," new ptarray required");
600  /* We didn't add the previous point, so this is a new segment.
601  * Make a new point array. */
602  dp = ptarray_construct_empty(hasz, hasm, 32);
603 
604  /* We're transiting into the range so add an interpolated
605  * point at the range boundary.
606  * If we're on a boundary and crossing from the far side,
607  * we also need an interpolated point. */
608  if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
609  ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
610  ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
611  ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
612  {
613  double interpolation_value;
614  (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
615  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
617  LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
618  }
619  }
620  /* Add the current vertex to the point array. */
622  if ( ordinate_value_p == from || ordinate_value_p == to )
623  {
624  added_last_point = 2; /* Added on boundary. */
625  }
626  else
627  {
628  added_last_point = 1; /* Added inside range. */
629  }
630  }
631  /* Is this point inside the ordinate range? No. */
632  else
633  {
634  LWDEBUGF(4, " added_last_point (%d)", added_last_point);
635  if ( added_last_point == 1 )
636  {
637  /* We're transiting out of the range, so add an interpolated point
638  * to the point array at the range boundary. */
639  double interpolation_value;
640  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
641  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
643  LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
644  }
645  else if ( added_last_point == 2 )
646  {
647  /* We're out and the last point was on the boundary.
648  * If the last point was the near boundary, nothing to do.
649  * If it was the far boundary, we need an interpolated point. */
650  if ( from != to && (
651  (ordinate_value_q == from && ordinate_value_p > from) ||
652  (ordinate_value_q == to && ordinate_value_p < to) ) )
653  {
654  double interpolation_value;
655  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
656  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
658  LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
659  }
660  }
661  else if ( i && ordinate_value_q < from && ordinate_value_p > to )
662  {
663  /* We just hopped over the whole range, from bottom to top,
664  * so we need to add *two* interpolated points! */
665  dp = ptarray_construct(hasz, hasm, 2);
666  /* Interpolate lower point. */
667  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
668  ptarray_set_point4d(dp, 0, r);
669  /* Interpolate upper point. */
670  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
671  ptarray_set_point4d(dp, 1, r);
672  }
673  else if ( i && ordinate_value_q > to && ordinate_value_p < from )
674  {
675  /* We just hopped over the whole range, from top to bottom,
676  * so we need to add *two* interpolated points! */
677  dp = ptarray_construct(hasz, hasm, 2);
678  /* Interpolate upper point. */
679  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
680  ptarray_set_point4d(dp, 0, r);
681  /* Interpolate lower point. */
682  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
683  ptarray_set_point4d(dp, 1, r);
684  }
685  /* We have an extant point-array, save it out to a multi-line. */
686  if ( dp )
687  {
688  LWDEBUG(4, "saving pointarray to multi-line (1)");
689 
690  /* Only one point, so we have to make an lwpoint to hold this
691  * and set the overall output type to a generic collection. */
692  if ( dp->npoints == 1 )
693  {
694  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
695  lwgeom_out->type = COLLECTIONTYPE;
696  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
697 
698  }
699  else
700  {
701  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
702  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
703  }
704 
705  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
706  dp = NULL;
707  }
708  added_last_point = 0;
709 
710  }
711  }
712 
713  /* Still some points left to be saved out. */
714  if ( dp && dp->npoints > 0 )
715  {
716  LWDEBUG(4, "saving pointarray to multi-line (2)");
717  LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
718  LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
719 
720  if ( dp->npoints == 1 )
721  {
722  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
723  lwgeom_out->type = COLLECTIONTYPE;
724  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
725  }
726  else
727  {
728  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
729  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
730  }
731 
732  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
733  dp = NULL;
734  }
735 
736  lwfree(p);
737  lwfree(q);
738  lwfree(r);
739 
740  if ( lwgeom_out->bbox && lwgeom_out->ngeoms > 0 )
741  {
742  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
743  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
744  }
745 
746  return lwgeom_out;
747 
748 }
749 
751 lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
752 {
753  LWCOLLECTION *out_col;
754  LWCOLLECTION *out_offset;
755  int i;
756 
757  if ( ! lwin )
758  lwerror("lwgeom_clip_to_ordinate_range: null input geometry!");
759 
760  switch ( lwin->type )
761  {
762  case LINETYPE:
763  out_col = lwline_clip_to_ordinate_range((LWLINE*)lwin, ordinate, from, to);
764  break;
765  case MULTILINETYPE:
766  out_col = lwmline_clip_to_ordinate_range((LWMLINE*)lwin, ordinate, from, to);
767  break;
768  case MULTIPOINTTYPE:
769  out_col = lwmpoint_clip_to_ordinate_range((LWMPOINT*)lwin, ordinate, from, to);
770  break;
771  case POINTTYPE:
772  out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to);
773  break;
774  default:
775  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
776  return NULL;;
777  }
778 
779  /* Stop if result is NULL */
780  if ( out_col == NULL )
781  lwerror("lwgeom_clip_to_ordinate_range clipping routine returned NULL");
782 
783  /* Return if we aren't going to offset the result */
784  if ( FP_EQUALS(offset, 0.0) || lwgeom_is_empty(lwcollection_as_lwgeom(out_col)) )
785  return out_col;
786 
787  /* Construct a collection to hold our outputs. */
788  /* Things get ugly: GEOS offset drops Z's and M's so we have to drop ours */
789  out_offset = lwcollection_construct_empty(MULTILINETYPE, lwin->srid, 0, 0);
790 
791  /* Try and offset the linear portions of the return value */
792  for ( i = 0; i < out_col->ngeoms; i++ )
793  {
794  int type = out_col->geoms[i]->type;
795  if ( type == POINTTYPE )
796  {
797  lwnotice("lwgeom_clip_to_ordinate_range cannot offset a clipped point");
798  continue;
799  }
800  else if ( type == LINETYPE )
801  {
802  /* lwgeom_offsetcurve(line, offset, quadsegs, joinstyle (round), mitrelimit) */
803  LWGEOM *lwoff = lwgeom_offsetcurve(lwgeom_as_lwline(out_col->geoms[i]), offset, 8, 1, 5.0);
804  if ( ! lwoff )
805  {
806  lwerror("lwgeom_offsetcurve returned null");
807  }
808  lwcollection_add_lwgeom(out_offset, lwoff);
809  }
810  else
811  {
812  lwerror("lwgeom_clip_to_ordinate_range found an unexpected type (%s) in the offset routine",lwtype_name(type));
813  }
814  }
815 
816  return out_offset;
817 }
818 
820 lwgeom_locate_between(const LWGEOM *lwin, double from, double to, double offset)
821 {
822  if ( ! lwgeom_has_m(lwin) )
823  lwerror("Input geometry does not have a measure dimension");
824 
825  return lwgeom_clip_to_ordinate_range(lwin, 'M', from, to, offset);
826 }
827 
828 double
829 lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
830 {
831  POINT4D p, p_proj;
832  double ret = 0.0;
833 
834  if ( ! lwin )
835  lwerror("lwgeom_interpolate_point: null input geometry!");
836 
837  if ( ! lwgeom_has_m(lwin) )
838  lwerror("Input geometry does not have a measure dimension");
839 
840  if ( lwgeom_is_empty(lwin) || lwpoint_is_empty(lwpt) )
841  lwerror("Input geometry is empty");
842 
843  switch ( lwin->type )
844  {
845  case LINETYPE:
846  {
847  LWLINE *lwline = lwgeom_as_lwline(lwin);
848  lwpoint_getPoint4d_p(lwpt, &p);
849  ret = ptarray_locate_point(lwline->points, &p, NULL, &p_proj);
850  ret = p_proj.m;
851  break;
852  }
853  default:
854  lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
855  }
856  return ret;
857 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:501
double x
Definition: liblwgeom.h:308
#define LINETYPE
Definition: liblwgeom.h:61
double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
Find the measure value at the location on the line closest to the point.
LWLINE * lwline_measured_from_lwline(const LWLINE *lwline, double m_start, double m_end)
Add a measure dimension to a line, interpolating linearly from the start to the end value...
Definition: lwline.c:366
uint8_t type
Definition: liblwgeom.h:433
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...
Definition: ptarray.c:49
double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
Given a POINT4D and an ordinate number, return the value of the ordinate.
LWCOLLECTION * lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
Clip an input POINT between two values, on any ordinate input.
double m
Definition: liblwgeom.h:308
LWCOLLECTION * lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
Clip an input MULTIPOINT between two values, on any ordinate input.
char * r
Definition: cu_in_wkt.c:25
void lwfree(void *mem)
Definition: lwutil.c:190
static LWMPOINT * lwline_locate_along(const LWLINE *lwline, double m, double offset)
int npoints
Definition: liblwgeom.h:327
uint8_t type
Definition: liblwgeom.h:459
static int segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn)
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:425
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:57
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:315
#define MULTIPOINTTYPE
Definition: liblwgeom.h:63
LWMPOINT * lwmpoint_construct(int srid, const POINTARRAY *pa)
Definition: lwmpoint.c:39
LWGEOM * lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
void lwline_free(LWLINE *line)
Definition: lwline.c:63
GBOX * bbox
Definition: liblwgeom.h:461
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:778
LWGEOM * lwmline_as_lwgeom(const LWMLINE *obj)
Definition: lwgeom.c:209
int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
Given two points, a dimensionality, an ordinate, and an interpolation value generate a new point that...
int32_t srid
Definition: liblwgeom.h:377
#define FP_MIN(A, B)
LWCOLLECTION * lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
Take in a LINESTRING and return a MULTILINESTRING of those portions of the LINESTRING between the fro...
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:792
static POINTARRAY * ptarray_locate_along(const POINTARRAY *pa, double m, double offset)
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:542
int32_t srid
Definition: liblwgeom.h:355
int ngeoms
Definition: liblwgeom.h:437
LWCOLLECTION * lwgeom_locate_between(const LWGEOM *lwin, double from, double to, double offset)
Determine the segments along a measured line that fall within the m-range given.
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
uint8_t flags
Definition: liblwgeom.h:460
double ptarray_locate_point(const POINTARRAY *pa, const POINT4D *pt, double *dist, POINT4D *p_located)
Definition: ptarray.c:1267
#define FLAGS_SET_Z(flags, value)
Definition: liblwgeom.h:112
void lwmpoint_free(LWMPOINT *mpt)
Definition: lwmpoint.c:59
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:249
int ptarray_append_point(POINTARRAY *pa, const POINT4D *pt, int allow_duplicates)
Append a point to the end of an existing POINTARRAY If allow_duplicate is LW_TRUE, then a duplicate point will not be added.
Definition: ptarray.c:141
#define LW_FALSE
Definition: liblwgeom.h:52
LWPOINT * lwpoint_clone(const LWPOINT *lwgeom)
Definition: lwpoint.c:206
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
LWGEOM ** geoms
Definition: liblwgeom.h:465
void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
Given a point, ordinate number and value, set that ordinate on the point.
static LWMPOINT * lwmline_locate_along(const LWMLINE *lwmline, double m, double offset)
int lwpoint_getPoint4d_p(const LWPOINT *point, POINT4D *out)
Definition: lwpoint.c:42
LWPOINT ** geoms
Definition: liblwgeom.h:426
LWMPOINT * lwmpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwmpoint.c:26
LWCOLLECTION * lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
Clip an input MULTILINESTRING between two values, on any ordinate input.
double z
Definition: liblwgeom.h:308
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:89
int32_t srid
Definition: liblwgeom.h:366
LWGEOM * lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
Determine the location(s) along a measured line where m occurs and return as a multipoint.
LWLINE ** geoms
Definition: liblwgeom.h:439
double lwpoint_get_m(const LWPOINT *point)
Definition: lwpoint.c:80
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:30
static LWMPOINT * lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset)
LWMPOINT * lwmpoint_add_lwpoint(LWMPOINT *mobj, const LWPOINT *obj)
Definition: lwmpoint.c:32
void * lwrealloc(void *mem, size_t size)
Definition: lwutil.c:183
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
int lwpoint_is_empty(const LWPOINT *point)
Definition: lwpoint.c:258
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:555
LWCOLLECTION * lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
Given a geometry clip based on the from/to range of one of its ordinates (x, y, z, m).
uint8_t type
Definition: liblwgeom.h:352
#define FP_EQUALS(A, B)
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:23
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:96
void * lwalloc(size_t size)
Definition: lwutil.c:175
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1229
double y
Definition: liblwgeom.h:308
#define MULTILINETYPE
Definition: liblwgeom.h:64
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
int ngeoms
Definition: liblwgeom.h:424
uint8_t flags
Definition: liblwgeom.h:375
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:118
int32_t srid
Definition: liblwgeom.h:423
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:799
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:174
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:217
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
static LWMPOINT * lwpoint_locate_along(const LWPOINT *lwpoint, double m, double offset)
#define FP_MAX(A, B)
int32_t srid
Definition: liblwgeom.h:436
#define FLAGS_SET_M(flags, value)
Definition: liblwgeom.h:113
POINTARRAY * points
Definition: liblwgeom.h:378
LWGEOM * lwcollection_as_lwgeom(const LWCOLLECTION *obj)
Definition: lwgeom.c:219
LWGEOM * lwmpoint_as_lwgeom(const LWMPOINT *obj)
Definition: lwgeom.c:214