lwgeom_distance_spheroid()

double lwgeom_distance_spheroid	(	const LWGEOM *	lwgeom1,
		const LWGEOM *	lwgeom2,
		const SPHEROID *	spheroid,
		double	tolerance
	)

Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid.

A spheroid with major axis == minor axis will be treated as a sphere. Pass in a tolerance in spheroid units.

Calculate the geodetic distance from lwgeom1 to lwgeom2 on the spheroid.

Return immediately when the calulated distance drops below the tolerance (useful for dwithin calculations). Return a negative distance for incalculable cases.

Definition at line 2183 of file lwgeodetic.c.

References LWGEOM::bbox, distance(), gbox_init(), gbox_overlaps(), LWCOLLECTION::geoms, getPoint2d_cp(), LINETYPE, LW_FALSE, LW_TRUE, LWDEBUG, LWDEBUGF, lwerror(), lwgeom_calculate_gbox_geodetic(), lwgeom_is_empty(), lwpoly_covers_point2d(), lwtype_is_collection(), lwtype_name(), LWCOLLECTION::ngeoms, LWPOLY::nrings, LWPOINT::point, LWLINE::points, POINTTYPE, POLYGONTYPE, ptarray_distance_spheroid(), LWPOLY::rings, and LWGEOM::type.

Referenced by geography_centroid_from_mline(), geography_distance(), geography_distance_knn(), geography_distance_uncached(), geography_dwithin(), geography_dwithin_uncached(), geometry_distance_spheroid(), test_lwgeom_distance_sphere(), test_tree_circ_distance(), and test_tree_circ_distance_threshold().

{
        uint8_t type1, type2;
        int check_intersection = LW_FALSE;
        GBOX gbox1, gbox2;

        gbox_init(&gbox1);
        gbox_init(&gbox2);

        assert(lwgeom1);
        assert(lwgeom2);

        LWDEBUGF(4, "entered function, tolerance %.8g", tolerance);

        /* What's the distance to an empty geometry? We don't know.
           Return a negative number so the caller can catch this case. */
        if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
        {
                return -1.0;
        }

        type1 = lwgeom1->type;
        type2 = lwgeom2->type;

        /* Make sure we have boxes */
        if ( lwgeom1->bbox )
                gbox1 = *(lwgeom1->bbox);
        else
                lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);

        /* Make sure we have boxes */
        if ( lwgeom2->bbox )
                gbox2 = *(lwgeom2->bbox);
        else
                lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);

        /* If the boxes aren't disjoint, we have to check for edge intersections */
        if ( gbox_overlaps(&gbox1, &gbox2) )
                check_intersection = LW_TRUE;

        /* Point/line combinations can all be handled with simple point array iterations */
        if ( ( type1 == POINTTYPE || type1 == LINETYPE ) &&
             ( type2 == POINTTYPE || type2 == LINETYPE ) )
        {
                POINTARRAY *pa1, *pa2;

                if ( type1 == POINTTYPE )
                        pa1 = ((LWPOINT*)lwgeom1)->point;
                else
                        pa1 = ((LWLINE*)lwgeom1)->points;

                if ( type2 == POINTTYPE )
                        pa2 = ((LWPOINT*)lwgeom2)->point;
                else
                        pa2 = ((LWLINE*)lwgeom2)->points;

                return ptarray_distance_spheroid(pa1, pa2, spheroid, tolerance, check_intersection);
        }

        /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */
        if ( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) ||
             ( type2 == POLYGONTYPE && type1 == POINTTYPE ) )
        {
                const POINT2D *p;
                LWPOLY *lwpoly;
                LWPOINT *lwpt;
                double distance = FLT_MAX;
                int i;

                if ( type1 == POINTTYPE )
                {
                        lwpt = (LWPOINT*)lwgeom1;
                        lwpoly = (LWPOLY*)lwgeom2;
                }
                else
                {
                        lwpt = (LWPOINT*)lwgeom2;
                        lwpoly = (LWPOLY*)lwgeom1;
                }
                p = getPoint2d_cp(lwpt->point, 0);

                /* Point in polygon implies zero distance */
                if ( lwpoly_covers_point2d(lwpoly, p) )
                {
                        return 0.0;
                }

                /* Not inside, so what's the actual distance? */
                for ( i = 0; i < lwpoly->nrings; i++ )
                {
                        double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, spheroid, tolerance, check_intersection);
                        if ( ring_distance < distance )
                                distance = ring_distance;
                        if ( distance < tolerance )
                                return distance;
                }
                return distance;
        }

        /* Line/polygon case, if start point-in-poly, return zero, else return distance. */
        if ( ( type1 == POLYGONTYPE && type2 == LINETYPE ) ||
             ( type2 == POLYGONTYPE && type1 == LINETYPE ) )
        {
                const POINT2D *p;
                LWPOLY *lwpoly;
                LWLINE *lwline;
                double distance = FLT_MAX;
                int i;

                if ( type1 == LINETYPE )
                {
                        lwline = (LWLINE*)lwgeom1;
                        lwpoly = (LWPOLY*)lwgeom2;
                }
                else
                {
                        lwline = (LWLINE*)lwgeom2;
                        lwpoly = (LWPOLY*)lwgeom1;
                }
                p = getPoint2d_cp(lwline->points, 0);

                LWDEBUG(4, "checking if a point of line is in polygon");

                /* Point in polygon implies zero distance */
                if ( lwpoly_covers_point2d(lwpoly, p) )
                        return 0.0;

                LWDEBUG(4, "checking ring distances");

                /* Not contained, so what's the actual distance? */
                for ( i = 0; i < lwpoly->nrings; i++ )
                {
                        double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, spheroid, tolerance, check_intersection);
                        LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance);
                        if ( ring_distance < distance )
                                distance = ring_distance;
                        if ( distance < tolerance )
                                return distance;
                }
                LWDEBUGF(4, "all rings checked, returning distance = %.8g", distance);
                return distance;

        }

        /* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */
        if ( ( type1 == POLYGONTYPE && type2 == POLYGONTYPE ) ||
             ( type2 == POLYGONTYPE && type1 == POLYGONTYPE ) )
        {
                const POINT2D *p;
                LWPOLY *lwpoly1 = (LWPOLY*)lwgeom1;
                LWPOLY *lwpoly2 = (LWPOLY*)lwgeom2;
                double distance = FLT_MAX;
                int i, j;

                /* Point of 2 in polygon 1 implies zero distance */
                p = getPoint2d_cp(lwpoly1->rings[0], 0);
                if ( lwpoly_covers_point2d(lwpoly2, p) )
                        return 0.0;

                /* Point of 1 in polygon 2 implies zero distance */
                p = getPoint2d_cp(lwpoly2->rings[0], 0);
                if ( lwpoly_covers_point2d(lwpoly1, p) )
                        return 0.0;

                /* Not contained, so what's the actual distance? */
                for ( i = 0; i < lwpoly1->nrings; i++ )
                {
                        for ( j = 0; j < lwpoly2->nrings; j++ )
                        {
                                double ring_distance = ptarray_distance_spheroid(lwpoly1->rings[i], lwpoly2->rings[j], spheroid, tolerance, check_intersection);
                                if ( ring_distance < distance )
                                        distance = ring_distance;
                                if ( distance < tolerance )
                                        return distance;
                        }
                }
                return distance;
        }

        /* Recurse into collections */
        if ( lwtype_is_collection(type1) )
        {
                int i;
                double distance = FLT_MAX;
                LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;

                for ( i = 0; i < col->ngeoms; i++ )
                {
                        double geom_distance = lwgeom_distance_spheroid(col->geoms[i], lwgeom2, spheroid, tolerance);
                        if ( geom_distance < distance )
                                distance = geom_distance;
                        if ( distance < tolerance )
                                return distance;
                }
                return distance;
        }

        /* Recurse into collections */
        if ( lwtype_is_collection(type2) )
        {
                int i;
                double distance = FLT_MAX;
                LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;

                for ( i = 0; i < col->ngeoms; i++ )
                {
                        double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], spheroid, tolerance);
                        if ( geom_distance < distance )
                                distance = geom_distance;
                        if ( distance < tolerance )
                                return distance;
                }
                return distance;
        }


        lwerror("arguments include unsupported geometry type (%s, %s)", lwtype_name(type1), lwtype_name(type1));
        return -1.0;

}

Here is the call graph for this function:

Here is the caller graph for this function: