PostGIS  2.3.7dev-r@@SVN_REVISION@@
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 from/to range for the specified ordinate (XYZM)

Clip a line based on the from/to range of one of its ordinates.

Definition at line 551 of file lwlinearreferencing.c.

References LWCOLLECTION::bbox, COLLECTIONTYPE, LWLINE::flags, FLAGS_NDIMS, getPoint4d_p(), LW_FALSE, lwalloc(), lwcollection_add_lwgeom(), lwcollection_construct_empty(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_add_bbox(), lwgeom_drop_bbox(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_to_ewkt(), lwline_as_lwgeom(), lwline_construct(), lwpoint_as_lwgeom(), lwpoint_construct(), lwpoint_get_ordinate(), MULTILINETYPE, LWCOLLECTION::ngeoms, POINTARRAY::npoints, point_interpolate(), LWLINE::points, ptarray_append_point(), ptarray_construct(), ptarray_construct_empty(), ptarray_set_point4d(), r, LWLINE::srid, and LWCOLLECTION::type.

Referenced by lwgeom_clip_to_ordinate_range(), lwmline_clip_to_ordinate_range(), test_lwline_clip(), test_lwline_clip_big(), and test_lwmline_clip().

552 {
553 
554  POINTARRAY *pa_in = NULL;
555  LWCOLLECTION *lwgeom_out = NULL;
556  POINTARRAY *dp = NULL;
557  int i;
558  int added_last_point = 0;
559  POINT4D *p = NULL, *q = NULL, *r = NULL;
560  double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
561  char hasz = lwgeom_has_z(lwline_as_lwgeom(line));
562  char hasm = lwgeom_has_m(lwline_as_lwgeom(line));
563  char dims = FLAGS_NDIMS(line->flags);
564 
565  /* Null input, nothing we can do. */
566  if ( ! line )
567  {
568  lwerror("Null input geometry.");
569  return NULL;
570  }
571 
572  /* Ensure 'from' is less than 'to'. */
573  if ( to < from )
574  {
575  double t = from;
576  from = to;
577  to = t;
578  }
579 
580  LWDEBUGF(4, "from = %g, to = %g, ordinate = %c", from, to, ordinate);
581  LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
582 
583  /* Asking for an ordinate we don't have. Error. */
584  if ( (ordinate == 'Z' && ! hasz) || (ordinate == 'M' && ! hasm) )
585  {
586  lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
587  return NULL;
588  }
589 
590  /* Prepare our working point objects. */
591  p = lwalloc(sizeof(POINT4D));
592  q = lwalloc(sizeof(POINT4D));
593  r = lwalloc(sizeof(POINT4D));
594 
595  /* Construct a collection to hold our outputs. */
596  lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
597 
598  /* Get our input point array */
599  pa_in = line->points;
600 
601  for ( i = 0; i < pa_in->npoints; i++ )
602  {
603  LWDEBUGF(4, "Point #%d", i);
604  LWDEBUGF(4, "added_last_point %d", added_last_point);
605  if ( i > 0 )
606  {
607  *q = *p;
608  ordinate_value_q = ordinate_value_p;
609  }
610  getPoint4d_p(pa_in, i, p);
611  ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
612  LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
613  LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
614 
615  /* Is this point inside the ordinate range? Yes. */
616  if ( ordinate_value_p >= from && ordinate_value_p <= to )
617  {
618  LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
619 
620  if ( ! added_last_point )
621  {
622  LWDEBUG(4," new ptarray required");
623  /* We didn't add the previous point, so this is a new segment.
624  * Make a new point array. */
625  dp = ptarray_construct_empty(hasz, hasm, 32);
626 
627  /* We're transiting into the range so add an interpolated
628  * point at the range boundary.
629  * If we're on a boundary and crossing from the far side,
630  * we also need an interpolated point. */
631  if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
632  ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
633  ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
634  ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
635  {
636  double interpolation_value;
637  (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
638  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
640  LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
641  }
642  }
643  /* Add the current vertex to the point array. */
645  if ( ordinate_value_p == from || ordinate_value_p == to )
646  {
647  added_last_point = 2; /* Added on boundary. */
648  }
649  else
650  {
651  added_last_point = 1; /* Added inside range. */
652  }
653  }
654  /* Is this point inside the ordinate range? No. */
655  else
656  {
657  LWDEBUGF(4, " added_last_point (%d)", added_last_point);
658  if ( added_last_point == 1 )
659  {
660  /* We're transiting out of the range, so add an interpolated point
661  * to the point array at the range boundary. */
662  double interpolation_value;
663  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
664  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
666  LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
667  }
668  else if ( added_last_point == 2 )
669  {
670  /* We're out and the last point was on the boundary.
671  * If the last point was the near boundary, nothing to do.
672  * If it was the far boundary, we need an interpolated point. */
673  if ( from != to && (
674  (ordinate_value_q == from && ordinate_value_p > from) ||
675  (ordinate_value_q == to && ordinate_value_p < to) ) )
676  {
677  double interpolation_value;
678  (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
679  point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
681  LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
682  }
683  }
684  else if ( i && ordinate_value_q < from && ordinate_value_p > to )
685  {
686  /* We just hopped over the whole range, from bottom to top,
687  * so we need to add *two* interpolated points! */
688  dp = ptarray_construct(hasz, hasm, 2);
689  /* Interpolate lower point. */
690  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
691  ptarray_set_point4d(dp, 0, r);
692  /* Interpolate upper point. */
693  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
694  ptarray_set_point4d(dp, 1, r);
695  }
696  else if ( i && ordinate_value_q > to && ordinate_value_p < from )
697  {
698  /* We just hopped over the whole range, from top to bottom,
699  * so we need to add *two* interpolated points! */
700  dp = ptarray_construct(hasz, hasm, 2);
701  /* Interpolate upper point. */
702  point_interpolate(p, q, r, hasz, hasm, ordinate, to);
703  ptarray_set_point4d(dp, 0, r);
704  /* Interpolate lower point. */
705  point_interpolate(p, q, r, hasz, hasm, ordinate, from);
706  ptarray_set_point4d(dp, 1, r);
707  }
708  /* We have an extant point-array, save it out to a multi-line. */
709  if ( dp )
710  {
711  LWDEBUG(4, "saving pointarray to multi-line (1)");
712 
713  /* Only one point, so we have to make an lwpoint to hold this
714  * and set the overall output type to a generic collection. */
715  if ( dp->npoints == 1 )
716  {
717  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
718  lwgeom_out->type = COLLECTIONTYPE;
719  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
720 
721  }
722  else
723  {
724  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
725  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
726  }
727 
728  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
729  dp = NULL;
730  }
731  added_last_point = 0;
732 
733  }
734  }
735 
736  /* Still some points left to be saved out. */
737  if ( dp && dp->npoints > 0 )
738  {
739  LWDEBUG(4, "saving pointarray to multi-line (2)");
740  LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
741  LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
742 
743  if ( dp->npoints == 1 )
744  {
745  LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
746  lwgeom_out->type = COLLECTIONTYPE;
747  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
748  }
749  else
750  {
751  LWLINE *oline = lwline_construct(line->srid, NULL, dp);
752  lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
753  }
754 
755  /* Pointarray is now owned by lwgeom_out, so drop reference to it */
756  dp = NULL;
757  }
758 
759  lwfree(p);
760  lwfree(q);
761  lwfree(r);
762 
763  if ( lwgeom_out->bbox && lwgeom_out->ngeoms > 0 )
764  {
765  lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
766  lwgeom_add_bbox((LWGEOM*)lwgeom_out);
767  }
768 
769  return lwgeom_out;
770 
771 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
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:62
double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
Given a POINT4D and an ordinate number, return the value of the ordinate.
char * r
Definition: cu_in_wkt.c:24
void lwfree(void *mem)
Definition: lwutil.c:242
int npoints
Definition: liblwgeom.h:370
uint8_t type
Definition: liblwgeom.h:502
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:482
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
GBOX * bbox
Definition: liblwgeom.h:504
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
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:420
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:599
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:262
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_FALSE, then a duplicate point will not be added.
Definition: ptarray.c:156
#define LW_FALSE
Definition: liblwgeom.h:76
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:612
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:267
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:111
void * lwalloc(size_t size)
Definition: lwutil.c:227
#define MULTILINETYPE
Definition: liblwgeom.h:88
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
uint8_t flags
Definition: liblwgeom.h:418
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:151
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:856
LWCOLLECTION * lwcollection_add_lwgeom(LWCOLLECTION *col, const LWGEOM *geom)
Appends geom to the collection managed by col.
Definition: lwcollection.c:187
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:244
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
POINTARRAY * points
Definition: liblwgeom.h:421

Here is the call graph for this function:

Here is the caller graph for this function: