lwhomogenize.o \
lwalgorithm.o \
lwsegmentize.o \
+ lwlinearreferencing.o \
lwprint.o \
vsprintf.o \
g_box.o \
extern double lwpoint_get_z(const LWPOINT *point);
extern double lwpoint_get_m(const LWPOINT *point);
+/**
+* Return SRID number
+*/
+extern int32_t lwgeom_get_srid(const LWGEOM *geom);
+
/**
* Return #LW_TRUE if geometry has Z ordinates
*/
*/
LWGEOM* lwgeom_sharedpaths(const LWGEOM* geom1, const LWGEOM* geom2);
+/*
+ * An offset curve against the input line.
+ *
+ * @param lwline a lineal geometry
+ * @param size offset distance. Offset left if negative and right if positive
+ * @param quadsegs number of quadrature segments in curves (try 8)
+ * @param joinStyle (1 = round, 2 = mitre, 3 = bevel)
+ * @param mitreLimit (try 5.0)
+ * @return derived geometry
+ *
+ * Requires GEOS-3.2.0+
+ */
+LWLINE* lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit);
+
/*******************************************************************************
* PROJ4-dependent extra functions on LWGEOM
-/*
-** lwpoint_get_ordinate(point, ordinate) => double
-*/
-double lwpoint_get_ordinate(const POINT4D *p, int ordinate)
-{
- if ( ! p )
- {
- lwerror("Null input geometry.");
- return 0.0;
- }
-
- if ( ordinate > 3 || ordinate < 0 )
- {
- lwerror("Cannot extract ordinate %d.", ordinate);
- return 0.0;
- }
-
- if ( ordinate == 3 )
- return p->m;
- if ( ordinate == 2 )
- return p->z;
- if ( ordinate == 1 )
- return p->y;
-
- return p->x;
-
-}
-void lwpoint_set_ordinate(POINT4D *p, int ordinate, double value)
-{
- if ( ! p )
- {
- lwerror("Null input geometry.");
- return;
- }
-
- if ( ordinate > 3 || ordinate < 0 )
- {
- lwerror("Cannot extract ordinate %d.", ordinate);
- return;
- }
-
- LWDEBUGF(4, " setting ordinate %d to %g", ordinate, value);
-
- switch ( ordinate )
- {
- case 3:
- p->m = value;
- return;
- case 2:
- p->z = value;
- return;
- case 1:
- p->y = value;
- return;
- case 0:
- p->x = value;
- return;
- }
-}
-
-
-int lwpoint_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value)
-{
- double p1_value = lwpoint_get_ordinate(p1, ordinate);
- double p2_value = lwpoint_get_ordinate(p2, ordinate);
- double proportion;
- int i = 0;
-
- if ( ordinate < 0 || ordinate >= ndims )
- {
- lwerror("Ordinate (%d) is not within ndims (%d).", ordinate, ndims);
- return 0;
- }
-
- if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
- FP_MAX(p1_value, p2_value) < interpolation_value )
- {
- lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
- return 0;
- }
-
- proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
-
- for ( i = 0; i < ndims; i++ )
- {
- double newordinate = 0.0;
- p1_value = lwpoint_get_ordinate(p1, i);
- p2_value = lwpoint_get_ordinate(p2, i);
- newordinate = p1_value + proportion * (p2_value - p1_value);
- lwpoint_set_ordinate(p, i, newordinate);
- LWDEBUGF(4, " clip ordinate(%d) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", i, p1_value, p2_value, proportion, newordinate );
- }
-
- return 1;
-}
-
-LWCOLLECTION *lwmline_clip_to_ordinate_range(LWMLINE *mline, int ordinate, double from, double to)
-{
- LWCOLLECTION *lwgeom_out = NULL;
-
- if ( ! mline )
- {
- lwerror("Null input geometry.");
- return NULL;
- }
-
- if ( mline->ngeoms == 1)
- {
- lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
- }
- else
- {
- LWCOLLECTION *col;
- char hasz = FLAGS_GET_Z(mline->flags);
- char hasm = FLAGS_GET_M(mline->flags);
- int i, j;
- char homogeneous = 1;
- size_t geoms_size = 0;
- lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, mline->srid, hasz, hasm);
- FLAGS_SET_Z(lwgeom_out->flags, hasz);
- FLAGS_SET_M(lwgeom_out->flags, hasm);
- for ( i = 0; i < mline->ngeoms; i ++ )
- {
- col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
- if ( col )
- {
- /* Something was left after the clip. */
- if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
- {
- geoms_size += 16;
- if ( lwgeom_out->geoms )
- {
- lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
- }
- else
- {
- lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
- }
- }
- for ( j = 0; j < col->ngeoms; j++ )
- {
- lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
- lwgeom_out->ngeoms++;
- }
- if ( col->type != mline->type )
- {
- homogeneous = 0;
- }
- /* Shallow free the struct, leaving the geoms behind. */
- if ( col->bbox ) lwfree(col->bbox);
- lwfree(col->geoms);
- lwfree(col);
- }
- }
- lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
- lwgeom_add_bbox((LWGEOM*)lwgeom_out);
- if ( ! homogeneous )
- {
- lwgeom_out->type = COLLECTIONTYPE;
- }
- }
-
- if ( ! lwgeom_out || lwgeom_out->ngeoms == 0 ) /* Nothing left after clip. */
- {
- return NULL;
- }
-
- return lwgeom_out;
-
-}
-
-/*
-** lwline_clip_to_ordinate_range(line, ordinate, from, to) => lwmline
-**
-** Take in a LINESTRING and return a MULTILINESTRING of those portions of the
-** LINESTRING between the from/to range for the specified ordinate (XYZM)
-*/
-LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to)
-{
-
- POINTARRAY *pa_in = NULL;
- LWCOLLECTION *lwgeom_out = NULL;
- POINTARRAY *dp = NULL;
- int i, rv;
- int added_last_point = 0;
- POINT4D *p = NULL, *q = NULL, *r = NULL;
- double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
- char hasz = FLAGS_GET_Z(line->flags);
- char hasm = FLAGS_GET_M(line->flags);
- char dims = FLAGS_NDIMS(line->flags);
-
- /* Null input, nothing we can do. */
- if ( ! line )
- {
- lwerror("Null input geometry.");
- return NULL;
- }
-
- /* Ensure 'from' is less than 'to'. */
- if ( to < from )
- {
- double t = from;
- from = to;
- to = t;
- }
-
- LWDEBUGF(4, "from = %g, to = %g, ordinate = %d", from, to, ordinate);
- LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
-
- /* Asking for an ordinate we don't have. Error. */
- if ( ordinate >= dims )
- {
- lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
- return NULL;
- }
-
- /* Prepare our working point objects. */
- p = lwalloc(sizeof(POINT4D));
- q = lwalloc(sizeof(POINT4D));
- r = lwalloc(sizeof(POINT4D));
-
- /* Construct a collection to hold our outputs. */
- lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
-
- /* Get our input point array */
- pa_in = line->points;
-
- for ( i = 0; i < pa_in->npoints; i++ )
- {
- LWDEBUGF(4, "Point #%d", i);
- LWDEBUGF(4, "added_last_point %d", added_last_point);
- if ( i > 0 )
- {
- *q = *p;
- ordinate_value_q = ordinate_value_p;
- }
- rv = getPoint4d_p(pa_in, i, p);
- ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
- LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
- LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
-
- /* Is this point inside the ordinate range? Yes. */
- if ( ordinate_value_p >= from && ordinate_value_p <= to )
- {
- LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
-
- if ( ! added_last_point )
- {
- LWDEBUG(4," new ptarray required");
- /* We didn't add the previous point, so this is a new segment.
- * Make a new point array. */
- dp = ptarray_construct_empty(hasz, hasm, 32);
-
- /* We're transiting into the range so add an interpolated
- * point at the range boundary.
- * If we're on a boundary and crossing from the far side,
- * we also need an interpolated point. */
- if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
- ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
- ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
- ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
- {
- double interpolation_value;
- (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
- rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
- rv = ptarray_append_point(dp, r, LW_FALSE);
- LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
- }
- }
- /* Add the current vertex to the point array. */
- rv = ptarray_append_point(dp, p, LW_FALSE);
- if ( ordinate_value_p == from || ordinate_value_p == to )
- {
- added_last_point = 2; /* Added on boundary. */
- }
- else
- {
- added_last_point = 1; /* Added inside range. */
- }
- }
- /* Is this point inside the ordinate range? No. */
- else
- {
- LWDEBUGF(4, " added_last_point (%d)", added_last_point);
- if ( added_last_point == 1 )
- {
- /* We're transiting out of the range, so add an interpolated point
- * to the point array at the range boundary. */
- double interpolation_value;
- (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
- rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
- rv = ptarray_append_point(dp, r, LW_FALSE);
- LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
- }
- else if ( added_last_point == 2 )
- {
- /* We're out and the last point was on the boundary.
- * If the last point was the near boundary, nothing to do.
- * If it was the far boundary, we need an interpolated point. */
- if ( from != to && (
- (ordinate_value_q == from && ordinate_value_p > from) ||
- (ordinate_value_q == to && ordinate_value_p < to) ) )
- {
- double interpolation_value;
- (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
- rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
- rv = ptarray_append_point(dp, r, LW_FALSE);
- LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
- }
- }
- else if ( i && ordinate_value_q < from && ordinate_value_p > to )
- {
- /* We just hopped over the whole range, from bottom to top,
- * so we need to add *two* interpolated points! */
- dp = ptarray_construct(hasz, hasm, 2);
- /* Interpolate lower point. */
- rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
- ptarray_set_point4d(dp, 0, r);
- /* Interpolate upper point. */
- rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
- ptarray_set_point4d(dp, 1, r);
- }
- else if ( i && ordinate_value_q > to && ordinate_value_p < from )
- {
- /* We just hopped over the whole range, from top to bottom,
- * so we need to add *two* interpolated points! */
- dp = ptarray_construct(hasz, hasm, 2);
- /* Interpolate upper point. */
- rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
- ptarray_set_point4d(dp, 0, r);
- /* Interpolate lower point. */
- rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
- ptarray_set_point4d(dp, 1, r);
- }
- /* We have an extant point-array, save it out to a multi-line. */
- if ( dp )
- {
- LWDEBUG(4, "saving pointarray to multi-line (1)");
-
- /* Only one point, so we have to make an lwpoint to hold this
- * and set the overall output type to a generic collection. */
- if ( dp->npoints == 1 )
- {
- LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
- lwgeom_out->type = COLLECTIONTYPE;
- lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
-
- }
- else
- {
- LWLINE *oline = lwline_construct(line->srid, NULL, dp);
- lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
- }
-
- /* Pointarray is now owned by lwgeom_out, so drop reference to it */
- dp = NULL;
- }
- added_last_point = 0;
-
- }
- }
-
- /* Still some points left to be saved out. */
- if ( dp && dp->npoints > 0 )
- {
- LWDEBUG(4, "saving pointarray to multi-line (2)");
- LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
- LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
-
- if ( dp->npoints == 1 )
- {
- LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
- lwgeom_out->type = COLLECTIONTYPE;
- lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
- }
- else
- {
- LWLINE *oline = lwline_construct(line->srid, NULL, dp);
- lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
- }
-
- /* Pointarray is now owned by lwgeom_out, so drop reference to it */
- dp = NULL;
- }
-
- lwfree(p);
- lwfree(q);
- lwfree(r);
-
- if ( lwgeom_out->ngeoms > 0 )
- {
- lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
- lwgeom_add_bbox((LWGEOM*)lwgeom_out);
- return lwgeom_out;
- }
-
- return NULL;
-
-}
static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
}
}
+int32_t
+lwgeom_get_srid(const LWGEOM *geom)
+{
+ if ( ! geom ) return SRID_UNKNOWN;
+ return geom->srid;
+}
+
int
lwgeom_has_z(const LWGEOM *geom)
{
return ver;
}
-
-
LWGEOM *
lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
{
g3 = GEOSSharedPaths(g1,g2);
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+
if (g3 == NULL)
{
- GEOSGeom_destroy(g1);
- GEOSGeom_destroy(g2);
lwerror("GEOSSharedPaths: %s", lwgeom_geos_errmsg);
return NULL;
}
- GEOSGeom_destroy(g1);
- GEOSGeom_destroy(g2);
-
GEOSSetSRID(g3, srid);
out = GEOS2LWGEOM(g3, is3d);
+ GEOSGeom_destroy(g3);
if (out == NULL)
{
- GEOSGeom_destroy(g3);
lwerror("GEOS2LWGEOM threw an error");
return NULL;
}
- GEOSGeom_destroy(g3);
return out;
#endif /* POSTGIS_GEOS_VERSION >= 33 */
}
+
+LWLINE*
+lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
+{
+#if POSTGIS_GEOS_VERSION < 32
+ lwerror("lwgeom_offsetcurve: GEOS 3.2 or higher required");
+#else
+ GEOSGeometry *g1, *g3;
+ LWGEOM *lwgeom_result;
+
+ initGEOS(lwnotice, lwgeom_geos_error);
+
+ g1 = (GEOSGeometry *)LWGEOM2GEOS(lwline);
+ if ( ! g1 )
+ {
+ lwerror("lwgeom_offsetcurve: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
+ return NULL;
+ }
+
+#if POSTGIS_GEOS_VERSION < 33
+ /* Size is always positive for GEOSSingleSidedBuffer, and a flag determines left/right */
+ g3 = GEOSSingleSidedBuffer(g1, size < 0 ? -size : size,
+ quadsegs, joinStyle, mitreLimit,
+ size < 0 ? 0 : 1);
+#else
+ g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
+#endif
+ /* Don't need input geometry anymore */
+ GEOSGeom_destroy(g1);
+
+ if (g3 == NULL)
+ {
+ GEOSGeom_destroy(g1);
+ lwerror("GEOSOffsetCurve: %s", lwgeom_geos_errmsg);
+ return NULL;
+ }
+
+ LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
+
+ GEOSSetSRID(g3, lwgeom_get_srid(lwline));
+
+ lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwline));
+ GEOSGeom_destroy(g3);
+
+ if (lwgeom_result == NULL)
+ {
+ lwerror("lwgeom_offsetcurve: GEOS2LWGEOM returned null");
+ return NULL;
+ }
+
+ return lwgeom_result;
+
+#endif /* POSTGIS_GEOS_VERSION < 32 */
+}
\ No newline at end of file
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2011 Paul Ramsey
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "liblwgeom_internal.h"
+#include "lwgeom_log.h"
+
+
+/**
+* Given a POINT4D and an ordinate number, return
+* the value of the ordinate.
+* @param p input point
+* @param ordinate number (1=x, 2=y, 3=z, 4=m)
+* @return d value at that ordinate
+*/
+double lwpoint_get_ordinate(const POINT4D *p, int ordinate)
+{
+ if ( ! p )
+ {
+ lwerror("Null input geometry.");
+ return 0.0;
+ }
+
+ if ( ordinate > 3 || ordinate < 0 )
+ {
+ lwerror("Cannot extract ordinate %d.", ordinate);
+ return 0.0;
+ }
+
+ if ( ordinate == 3 )
+ return p->m;
+ if ( ordinate == 2 )
+ return p->z;
+ if ( ordinate == 1 )
+ return p->y;
+
+ return p->x;
+
+}
+
+/**
+* Given a point, ordinate number and value, set that ordinate on the
+* point.
+*/
+void lwpoint_set_ordinate(POINT4D *p, int ordinate, double value)
+{
+ if ( ! p )
+ {
+ lwerror("Null input geometry.");
+ return;
+ }
+
+ if ( ordinate > 3 || ordinate < 0 )
+ {
+ lwerror("Cannot extract ordinate %d.", ordinate);
+ return;
+ }
+
+ LWDEBUGF(4, " setting ordinate %d to %g", ordinate, value);
+
+ switch ( ordinate )
+ {
+ case 3:
+ p->m = value;
+ return;
+ case 2:
+ p->z = value;
+ return;
+ case 1:
+ p->y = value;
+ return;
+ case 0:
+ p->x = value;
+ return;
+ }
+}
+
+/**
+* Given two points, a dimensionality, an ordinate, and an interpolation value
+* generate a new point that is proportionally between the input points,
+* using the values in the provided dimension as the scaling factors.
+*/
+int lwpoint_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value)
+{
+ double p1_value = lwpoint_get_ordinate(p1, ordinate);
+ double p2_value = lwpoint_get_ordinate(p2, ordinate);
+ double proportion;
+ int i = 0;
+
+ if ( ordinate < 0 || ordinate >= ndims )
+ {
+ lwerror("Ordinate (%d) is not within ndims (%d).", ordinate, ndims);
+ return 0;
+ }
+
+ if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
+ FP_MAX(p1_value, p2_value) < interpolation_value )
+ {
+ lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
+ return 0;
+ }
+
+ proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
+
+ for ( i = 0; i < ndims; i++ )
+ {
+ double newordinate = 0.0;
+ p1_value = lwpoint_get_ordinate(p1, i);
+ p2_value = lwpoint_get_ordinate(p2, i);
+ newordinate = p1_value + proportion * (p2_value - p1_value);
+ lwpoint_set_ordinate(p, i, newordinate);
+ LWDEBUGF(4, " clip ordinate(%d) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", i, p1_value, p2_value, proportion, newordinate );
+ }
+
+ return 1;
+}
+
+/**
+* Clip an input MULTILINESTRING between two values, on any ordinate input.
+*/
+LWCOLLECTION *lwmline_clip_to_ordinate_range(LWMLINE *mline, int ordinate, double from, double to)
+{
+ LWCOLLECTION *lwgeom_out = NULL;
+
+ if ( ! mline )
+ {
+ lwerror("Null input geometry.");
+ return NULL;
+ }
+
+ if ( mline->ngeoms == 1)
+ {
+ lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
+ }
+ else
+ {
+ LWCOLLECTION *col;
+ char hasz = FLAGS_GET_Z(mline->flags);
+ char hasm = FLAGS_GET_M(mline->flags);
+ int i, j;
+ char homogeneous = 1;
+ size_t geoms_size = 0;
+ lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, mline->srid, hasz, hasm);
+ FLAGS_SET_Z(lwgeom_out->flags, hasz);
+ FLAGS_SET_M(lwgeom_out->flags, hasm);
+ for ( i = 0; i < mline->ngeoms; i ++ )
+ {
+ col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
+ if ( col )
+ {
+ /* Something was left after the clip. */
+ if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
+ {
+ geoms_size += 16;
+ if ( lwgeom_out->geoms )
+ {
+ lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
+ }
+ else
+ {
+ lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
+ }
+ }
+ for ( j = 0; j < col->ngeoms; j++ )
+ {
+ lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
+ lwgeom_out->ngeoms++;
+ }
+ if ( col->type != mline->type )
+ {
+ homogeneous = 0;
+ }
+ /* Shallow free the struct, leaving the geoms behind. */
+ if ( col->bbox ) lwfree(col->bbox);
+ lwfree(col->geoms);
+ lwfree(col);
+ }
+ }
+ lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
+ lwgeom_add_bbox((LWGEOM*)lwgeom_out);
+ if ( ! homogeneous )
+ {
+ lwgeom_out->type = COLLECTIONTYPE;
+ }
+ }
+
+ if ( ! lwgeom_out || lwgeom_out->ngeoms == 0 ) /* Nothing left after clip. */
+ {
+ return NULL;
+ }
+
+ return lwgeom_out;
+
+}
+
+
+/**
+* Take in a LINESTRING and return a MULTILINESTRING of those portions of the
+* LINESTRING between the from/to range for the specified ordinate (XYZM)
+*/
+LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to)
+{
+
+ POINTARRAY *pa_in = NULL;
+ LWCOLLECTION *lwgeom_out = NULL;
+ POINTARRAY *dp = NULL;
+ int i, rv;
+ int added_last_point = 0;
+ POINT4D *p = NULL, *q = NULL, *r = NULL;
+ double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
+ char hasz = FLAGS_GET_Z(line->flags);
+ char hasm = FLAGS_GET_M(line->flags);
+ char dims = FLAGS_NDIMS(line->flags);
+
+ /* Null input, nothing we can do. */
+ if ( ! line )
+ {
+ lwerror("Null input geometry.");
+ return NULL;
+ }
+
+ /* Ensure 'from' is less than 'to'. */
+ if ( to < from )
+ {
+ double t = from;
+ from = to;
+ to = t;
+ }
+
+ LWDEBUGF(4, "from = %g, to = %g, ordinate = %d", from, to, ordinate);
+ LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line));
+
+ /* Asking for an ordinate we don't have. Error. */
+ if ( ordinate >= dims )
+ {
+ lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
+ return NULL;
+ }
+
+ /* Prepare our working point objects. */
+ p = lwalloc(sizeof(POINT4D));
+ q = lwalloc(sizeof(POINT4D));
+ r = lwalloc(sizeof(POINT4D));
+
+ /* Construct a collection to hold our outputs. */
+ lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
+
+ /* Get our input point array */
+ pa_in = line->points;
+
+ for ( i = 0; i < pa_in->npoints; i++ )
+ {
+ LWDEBUGF(4, "Point #%d", i);
+ LWDEBUGF(4, "added_last_point %d", added_last_point);
+ if ( i > 0 )
+ {
+ *q = *p;
+ ordinate_value_q = ordinate_value_p;
+ }
+ rv = getPoint4d_p(pa_in, i, p);
+ ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
+ LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
+ LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
+
+ /* Is this point inside the ordinate range? Yes. */
+ if ( ordinate_value_p >= from && ordinate_value_p <= to )
+ {
+ LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
+
+ if ( ! added_last_point )
+ {
+ LWDEBUG(4," new ptarray required");
+ /* We didn't add the previous point, so this is a new segment.
+ * Make a new point array. */
+ dp = ptarray_construct_empty(hasz, hasm, 32);
+
+ /* We're transiting into the range so add an interpolated
+ * point at the range boundary.
+ * If we're on a boundary and crossing from the far side,
+ * we also need an interpolated point. */
+ if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
+ ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
+ ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
+ ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
+ {
+ double interpolation_value;
+ (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
+ rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
+ rv = ptarray_append_point(dp, r, LW_FALSE);
+ LWDEBUGF(4, "[0] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
+ }
+ }
+ /* Add the current vertex to the point array. */
+ rv = ptarray_append_point(dp, p, LW_FALSE);
+ if ( ordinate_value_p == from || ordinate_value_p == to )
+ {
+ added_last_point = 2; /* Added on boundary. */
+ }
+ else
+ {
+ added_last_point = 1; /* Added inside range. */
+ }
+ }
+ /* Is this point inside the ordinate range? No. */
+ else
+ {
+ LWDEBUGF(4, " added_last_point (%d)", added_last_point);
+ if ( added_last_point == 1 )
+ {
+ /* We're transiting out of the range, so add an interpolated point
+ * to the point array at the range boundary. */
+ double interpolation_value;
+ (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
+ rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
+ rv = ptarray_append_point(dp, r, LW_FALSE);
+ LWDEBUGF(4, " [1] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
+ }
+ else if ( added_last_point == 2 )
+ {
+ /* We're out and the last point was on the boundary.
+ * If the last point was the near boundary, nothing to do.
+ * If it was the far boundary, we need an interpolated point. */
+ if ( from != to && (
+ (ordinate_value_q == from && ordinate_value_p > from) ||
+ (ordinate_value_q == to && ordinate_value_p < to) ) )
+ {
+ double interpolation_value;
+ (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
+ rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
+ rv = ptarray_append_point(dp, r, LW_FALSE);
+ LWDEBUGF(4, " [2] interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
+ }
+ }
+ else if ( i && ordinate_value_q < from && ordinate_value_p > to )
+ {
+ /* We just hopped over the whole range, from bottom to top,
+ * so we need to add *two* interpolated points! */
+ dp = ptarray_construct(hasz, hasm, 2);
+ /* Interpolate lower point. */
+ rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
+ ptarray_set_point4d(dp, 0, r);
+ /* Interpolate upper point. */
+ rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
+ ptarray_set_point4d(dp, 1, r);
+ }
+ else if ( i && ordinate_value_q > to && ordinate_value_p < from )
+ {
+ /* We just hopped over the whole range, from top to bottom,
+ * so we need to add *two* interpolated points! */
+ dp = ptarray_construct(hasz, hasm, 2);
+ /* Interpolate upper point. */
+ rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
+ ptarray_set_point4d(dp, 0, r);
+ /* Interpolate lower point. */
+ rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
+ ptarray_set_point4d(dp, 1, r);
+ }
+ /* We have an extant point-array, save it out to a multi-line. */
+ if ( dp )
+ {
+ LWDEBUG(4, "saving pointarray to multi-line (1)");
+
+ /* Only one point, so we have to make an lwpoint to hold this
+ * and set the overall output type to a generic collection. */
+ if ( dp->npoints == 1 )
+ {
+ LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
+ lwgeom_out->type = COLLECTIONTYPE;
+ lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
+
+ }
+ else
+ {
+ LWLINE *oline = lwline_construct(line->srid, NULL, dp);
+ lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
+ }
+
+ /* Pointarray is now owned by lwgeom_out, so drop reference to it */
+ dp = NULL;
+ }
+ added_last_point = 0;
+
+ }
+ }
+
+ /* Still some points left to be saved out. */
+ if ( dp && dp->npoints > 0 )
+ {
+ LWDEBUG(4, "saving pointarray to multi-line (2)");
+ LWDEBUGF(4, "dp->npoints == %d", dp->npoints);
+ LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
+
+ if ( dp->npoints == 1 )
+ {
+ LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
+ lwgeom_out->type = COLLECTIONTYPE;
+ lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
+ }
+ else
+ {
+ LWLINE *oline = lwline_construct(line->srid, NULL, dp);
+ lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
+ }
+
+ /* Pointarray is now owned by lwgeom_out, so drop reference to it */
+ dp = NULL;
+ }
+
+ lwfree(p);
+ lwfree(q);
+ lwfree(r);
+
+ if ( lwgeom_out->ngeoms > 0 )
+ {
+ lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
+ lwgeom_add_bbox((LWGEOM*)lwgeom_out);
+ return lwgeom_out;
+ }
+
+ return NULL;
+
+}
\ No newline at end of file
Datum isvalidreason(PG_FUNCTION_ARGS);
Datum isvaliddetail(PG_FUNCTION_ARGS);
Datum buffer(PG_FUNCTION_ARGS);
-Datum offsetcurve(PG_FUNCTION_ARGS);
Datum intersection(PG_FUNCTION_ARGS);
Datum convexhull(PG_FUNCTION_ARGS);
Datum topologypreservesimplify(PG_FUNCTION_ARGS);
Datum hausdorffdistance(PG_FUNCTION_ARGS)
{
#if POSTGIS_GEOS_VERSION < 32
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'ST_HausdorffDistance' function (3.2.0+ required)",
POSTGIS_GEOS_VERSION);
Datum hausdorffdistancedensify(PG_FUNCTION_ARGS)
{
#if POSTGIS_GEOS_VERSION < 32
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'ST_HausdorffDistance' function (3.2.0+ required)",
POSTGIS_GEOS_VERSION);
{
#if POSTGIS_GEOS_VERSION < 33
PG_RETURN_NULL();
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'GEOSUnaryUnion' function (3.3.0+ required)",
POSTGIS_GEOS_VERSION);
endCapStyle != DEFAULT_ENDCAP_STYLE ||
joinStyle != DEFAULT_JOIN_STYLE )
{
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"specifying a mitre limit != %d or styles different "
"from 'round' (needs 3.2 or higher)",
PG_RETURN_POINTER(result);
}
-PG_FUNCTION_INFO_V1(offsetcurve);
-Datum offsetcurve(PG_FUNCTION_ARGS)
+/*
+* Compute at offset curve to a line
+*/
+Datum ST_OffsetCurve(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_OffsetCurve);
+Datum ST_OffsetCurve(PG_FUNCTION_ARGS)
{
-#if POSTGIS_GEOS_VERSION >= 32
- GSERIALIZED *geom1;
- double size;
- GEOSGeometry *g1, *g3;
- GSERIALIZED *result;
+#if POSTGIS_GEOS_VERSION < 32
+ lwerror("The GEOS version this PostGIS binary "
+ "was compiled against (%d) doesn't support "
+ "ST_OffsetCurve function "
+ "(needs 3.2 or higher)",
+ POSTGIS_GEOS_VERSION);
+ PG_RETURN_NULL(); /* never get here */
+#else
+
+ GSERIALIZED *gser_input;
+ GSERIALIZED *gser_result;
+ LWGEOM *lwgeom_input;
+ LWGEOM *lwgeom_result;
+ double size;
int quadsegs = 8; /* the default */
int nargs;
JOIN_MITRE = 2,
JOIN_BEVEL = 3
};
+
static const double DEFAULT_MITRE_LIMIT = 5.0;
static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
-
double mitreLimit = DEFAULT_MITRE_LIMIT;
int joinStyle = DEFAULT_JOIN_STYLE;
- char *param;
- char *params = NULL;
-
-
- geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ char *param = NULL;
+ char *paramstr = NULL;
+
+ /* Read SQL arguments */
+ nargs = PG_NARGS();
+ gser_input = (GSERIALIZED*) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
size = PG_GETARG_FLOAT8(1);
- /*
- * For distance = 0 we just return the input.
- * Note that due to a bug, GEOS 3.3.0 would return EMPTY.
- * See http://trac.osgeo.org/geos/ticket/454
- */
- if ( size == 0 ) {
- PG_RETURN_POINTER(geom1);
+ /* Check for a useable type */
+ if ( gserialized_get_type(gser_input) != LINETYPE )
+ {
+ lwerror("ST_OffsetCurve only works with LineStrings");
+ PG_RETURN_NULL();
}
+ /*
+ * For distance == 0, just return the input.
+ * Note that due to a bug, GEOS 3.3.0 would return EMPTY.
+ * See http://trac.osgeo.org/geos/ticket/454
+ */
+ if ( size == 0 )
+ PG_RETURN_POINTER(gser_input);
+ /* Read the lwgeom, check for errors */
+ lwgeom_input = lwgeom_from_gserialized(gser_input);
+ if ( ! lwgeom_input )
+ lwerror("ST_OffsetCurve: lwgeom_from_gserialized returned NULL");
+
+ /* For empty inputs, just echo them back */
+ if ( lwgeom_is_empty(lwgeom_input) )
+ PG_RETURN_POINTER(gser_input);
- nargs = PG_NARGS();
-
- initGEOS(lwnotice, lwnotice);
- initGEOS(lwnotice, lwgeom_geos_error);
-
- g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
- if ( ! g1 ) {
- lwerror("Geometry could not be converted to GEOS: %s",
- lwgeom_geos_errmsg);
- PG_RETURN_NULL();
- }
-
- // options arg (optional)
- if (nargs > 2)
+ /* Process the optional arguments */
+ if ( nargs > 2 )
{
- /* We strdup `cause we're going to modify it */
- {
- text *wkttext = PG_GETARG_TEXT_P(2);
- params = text2cstring(wkttext);
- }
- /*params = pstrdup(PG_GETARG_CSTRING(2)); */
+ text *wkttext = PG_GETARG_TEXT_P(2);
+ paramstr = text2cstring(wkttext);
- POSTGIS_DEBUGF(3, "Params: %s", params);
+ POSTGIS_DEBUGF(3, "paramstr: %s", paramstr);
- for (param=params; ; param=NULL)
+ for ( param=paramstr; ; param=NULL )
{
char *key, *val;
param = strtok(param, " ");
val = strchr(key, '=');
if ( val == NULL || *(val+1) == '\0' )
{
- lwerror("Missing value for buffer "
- "parameter %s", key);
+ lwerror("ST_OffsetCurve: Missing value for buffer parameter %s", key);
break;
}
*val = '\0';
{
joinStyle = JOIN_ROUND;
}
- else if ( !strcmp(val, "mitre") ||
- !strcmp(val, "miter") )
+ else if ( !(strcmp(val, "mitre") && strcmp(val, "miter")) )
{
joinStyle = JOIN_MITRE;
}
- else if ( !strcmp(val, "bevel") )
+ else if ( ! strcmp(val, "bevel") )
{
joinStyle = JOIN_BEVEL;
}
else
{
- lwerror("Invalid buffer end cap "
- "style: %s (accept: "
- "'round', 'mitre', 'miter' "
- " or 'bevel'"
- ")", val);
+ lwerror("Invalid buffer end cap style: %s (accept: "
+ "'round', 'mitre', 'miter' or 'bevel')", val);
break;
}
}
else
{
lwerror("Invalid buffer parameter: %s (accept: "
- "'join', 'mitre_limit', "
- "'miter_limit and "
+ "'join', 'mitre_limit', 'miter_limit and "
"'quad_segs')", key);
break;
}
}
-
- pfree(params); /* was pstrduped */
-
- POSTGIS_DEBUGF(3, "joinStyle:%d mitreLimit:%g",
- joinStyle, mitreLimit);
-
+ POSTGIS_DEBUGF(3, "joinStyle:%d mitreLimit:%g", joinStyle, mitreLimit);
+ pfree(paramstr); /* alloc'ed in text2cstring */
}
+
+ lwgeom_result = lwline_as_lwgeom(lwgeom_offsetcurve(lwgeom_as_lwline(lwgeom_input), size, quadsegs, joinStyle, mitreLimit));
+
+ if (lwgeom_result == NULL)
+ lwerror("ST_OffsetCurve: lwgeom_offsetcurve returned NULL");
-#if POSTGIS_GEOS_VERSION < 33
- g3 = GEOSSingleSidedBuffer(g1, size < 0 ? -size : size,
- quadsegs, joinStyle, mitreLimit,
- size < 0 ? 0 : 1);
-#else
- g3 = GEOSOffsetCurve(g1, size, quadsegs, joinStyle, mitreLimit);
-#endif
-
- if (g3 == NULL)
- {
- lwerror("GEOSOffsetCurve: %s", lwgeom_geos_errmsg);
- GEOSGeom_destroy(g1);
- PG_RETURN_NULL(); /* never get here */
- }
-
- POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
-
- GEOSSetSRID(g3, gserialized_get_srid(geom1));
-
- result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
-
- if (result == NULL)
- {
- GEOSGeom_destroy(g1);
- GEOSGeom_destroy(g3);
- lwerror("ST_OffsetCurve() threw an error (result postgis geometry formation)!");
- PG_RETURN_NULL(); /* never get here */
- }
- GEOSGeom_destroy(g1);
- GEOSGeom_destroy(g3);
-
-
- /* compressType(result); */
-
- PG_FREE_IF_COPY(geom1, 0);
+ gser_result = gserialized_from_lwgeom(lwgeom_result, 0, 0);
+ lwgeom_free(lwgeom_input);
+ lwgeom_free(lwgeom_result);
+ PG_RETURN_POINTER(gser_result);
- PG_RETURN_POINTER(result);
-#else /* POSTGIS_GEOS_VERSION < 32 */
- lwerror("The GEOS version this postgis binary "
- "was compiled against (%d) doesn't support "
- "ST_OffsetCurve function "
- "(needs 3.2 or higher)",
- POSTGIS_GEOS_VERSION);
- PG_RETURN_NULL(); /* never get here */
#endif /* POSTGIS_GEOS_VERSION < 32 */
}
Datum isvaliddetail(PG_FUNCTION_ARGS)
{
#if POSTGIS_GEOS_VERSION < 33
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'isValidDetail' function (3.3.0+ required)",
POSTGIS_GEOS_VERSION);
#if POSTGIS_GEOS_VERSION >= 33
bnr = PG_GETARG_INT32(2);
#else
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"specifying a boundary node rule with ST_Relate"
" (3.3.0+ required)",
Datum ST_Snap(PG_FUNCTION_ARGS)
{
#if POSTGIS_GEOS_VERSION < 33
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'ST_Snap' function (3.3.0+ required)",
POSTGIS_GEOS_VERSION);
Datum ST_SharedPaths(PG_FUNCTION_ARGS)
{
#if POSTGIS_GEOS_VERSION < 33
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'ST_SharedPaths' function (3.3.0+ required)",
POSTGIS_GEOS_VERSION);
Datum ST_Node(PG_FUNCTION_ARGS)
{
#if POSTGIS_GEOS_VERSION < 33
- lwerror("The GEOS version this postgis binary "
+ lwerror("The GEOS version this PostGIS binary "
"was compiled against (%d) doesn't support "
"'ST_Node' function (3.3.0+ required)",
POSTGIS_GEOS_VERSION);
poObj = json_tokener_parse_ex(jstok, geojson, -1);
if( jstok->err != json_tokener_success)
{
- char *err;
- err += sprintf(err, "%s (at offset %d)", json_tokener_errors[jstok->err], jstok->char_offset);
+ char err[256];
+ snprintf(err, 256, "%s (at offset %d)", json_tokener_errors[jstok->err], jstok->char_offset);
json_tokener_free(jstok);
geojson_lwerror(err, 1);
}
-- Availability: 2.0.0 - requires GEOS-3.2 or higher\r
CREATE OR REPLACE FUNCTION ST_OffsetCurve(line geometry, distance float8, params text DEFAULT '')\r
RETURNS geometry\r
- AS 'MODULE_PATHNAME','offsetcurve'\r
+ AS 'MODULE_PATHNAME','ST_OffsetCurve'\r
LANGUAGE 'C' IMMUTABLE STRICT\r
COST 100;\r
\r
-ERROR: GEOSOffsetCurve: IllegalArgumentException: BufferBuilder::bufferLineSingleSided only accept linestrings
+ERROR: ST_OffsetCurve only works with LineStrings
t0|SRID=42;LINESTRING(0 0,10 0)
t1|SRID=42;LINESTRING(0 10,10 10)
t2|SRID=42;LINESTRING(10 -10,0 -10)
projects / postgis / commitdiff