geography_interpolate_points()

LWGEOM* geography_interpolate_points	(	const LWLINE *	line,
		double	length_fraction,
		const SPHEROID *	s,
		char	repeat
	)

Interpolate a point along a geographic line.

Definition at line 256 of file lwgeodetic_measures.c.

{
        POINT4D pt;
        uint32_t i;
        uint32_t points_to_interpolate;
        uint32_t points_found = 0;
        double length;
        double length_fraction_increment = length_fraction;
        double length_fraction_consumed = 0;
        char has_z = (char) lwgeom_has_z(lwline_as_lwgeom(line));
        char has_m = (char) lwgeom_has_m(lwline_as_lwgeom(line));
        const POINTARRAY* ipa = line->points;
        POINTARRAY *opa;
        POINT4D p1, p2;
        POINT3D q1, q2;
        LWGEOM *lwresult;
        GEOGRAPHIC_POINT g1, g2;
        uint32_t srid = line->srid;
 
        /* Empty.InterpolatePoint == Point Empty */
        if ( lwline_is_empty(line) )
        {
                return lwgeom_clone_deep(lwline_as_lwgeom(line));
        }
 
        /*
        * If distance is one of the two extremes, return the point on that
        * end rather than doing any computations
        */
        if ( length_fraction == 0.0 || length_fraction == 1.0 )
        {
                if ( length_fraction == 0.0 )
                        getPoint4d_p(ipa, 0, &pt);
                else
                        getPoint4d_p(ipa, ipa->npoints-1, &pt);
 
                opa = ptarray_construct(has_z, has_m, 1);
                ptarray_set_point4d(opa, 0, &pt);
                return lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
        }
 
        /* Interpolate points along the line */
        length = ptarray_length_spheroid(ipa, s);
        points_to_interpolate = repeat ? (uint32_t) floor(1 / length_fraction) : 1;
        opa = ptarray_construct(has_z, has_m, points_to_interpolate);
 
        getPoint4d_p(ipa, 0, &p1);
        geographic_point_init(p1.x, p1.y, &g1);
        for ( i = 0; i < ipa->npoints - 1 && points_found < points_to_interpolate; i++ )
        {
                double segment_length_frac;
                getPoint4d_p(ipa, i+1, &p2);
                geographic_point_init(p2.x, p2.y, &g2);
 
                /* Special sphere case */
                if ( s->a == s->b )
                        segment_length_frac = s->radius * sphere_distance(&g1, &g2) / length;
                /* Spheroid case */
                else
                        segment_length_frac = spheroid_distance(&g1, &g2, s) / length;
 
                /* If our target distance is before the total length we've seen
                * so far. create a new point some distance down the current
                * segment.
                */
                while ( length_fraction < length_fraction_consumed + segment_length_frac && points_found < points_to_interpolate )
                {
                        geog2cart(&g1, &q1);
                        geog2cart(&g2, &q2);
                        double segment_fraction = (length_fraction - length_fraction_consumed) / segment_length_frac;
                        interpolate_point4d_spheroid(&p1, &p2, &pt, s, segment_fraction);
                        ptarray_set_point4d(opa, points_found++, &pt);
                        length_fraction += length_fraction_increment;
                }
 
                length_fraction_consumed += segment_length_frac;
 
                p1 = p2;
                g1 = g2;
        }
 
        /* Return the last point on the line. This shouldn't happen, but
        * could if there's some floating point rounding errors. */
        if (points_found < points_to_interpolate)
        {
                getPoint4d_p(ipa, ipa->npoints - 1, &pt);
                ptarray_set_point4d(opa, points_found, &pt);
        }
 
        if (opa->npoints <= 1)
        {
                lwresult = lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
        } else {
                lwresult = lwmpoint_as_lwgeom(lwmpoint_construct(srid, opa));
        }
 
        return lwresult;
}

References geog2cart(), geographic_point_init(), getPoint4d_p(), interpolate_point4d_spheroid(), lwgeom_clone_deep(), lwgeom_has_m(), lwgeom_has_z(), lwline_as_lwgeom(), lwline_is_empty(), lwmpoint_as_lwgeom(), lwmpoint_construct(), lwpoint_as_lwgeom(), lwpoint_construct(), POINTARRAY::npoints, LWLINE::points, ptarray_construct(), ptarray_length_spheroid(), ptarray_set_point4d(), s, sphere_distance(), spheroid_distance(), LWLINE::srid, POINT4D::x, and POINT4D::y.

Referenced by geography_line_interpolate_point().

Here is the call graph for this function:

Here is the caller graph for this function: