PostGIS  2.5.7dev-r@@SVN_REVISION@@

◆ lwline_clip_to_ordinate_range()

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 547 of file lwlinearreferencing.c.

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

References LWCOLLECTION::bbox, COLLECTIONTYPE, LWLINE::flags, FLAGS_NDIMS, getPoint4d_p(), LW_FALSE, lwalloc(), lwcollection_add_lwgeom(), lwcollection_construct_empty(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_refresh_bbox(), 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().

Here is the call graph for this function:
Here is the caller graph for this function: