PostGIS  2.1.10dev-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 528 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().

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 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:501
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.
char * r
Definition: cu_in_wkt.c:25
void lwfree(void *mem)
Definition: lwutil.c:190
int npoints
Definition: liblwgeom.h:327
uint8_t type
Definition: liblwgeom.h:459
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
GBOX * bbox
Definition: liblwgeom.h:461
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
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
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:792
void lwgeom_drop_bbox(LWGEOM *lwgeom)
Call this function to drop BBOX and SRID from LWGEOM.
Definition: lwgeom.c:542
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
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
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:29
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:555
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:254
LWPOINT * lwpoint_construct(int srid, GBOX *bbox, POINTARRAY *point)
Definition: lwpoint.c:96
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
uint8_t flags
Definition: liblwgeom.h:375
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:118
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
POINTARRAY * points
Definition: liblwgeom.h:378

Here is the call graph for this function:

Here is the caller graph for this function: