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 calculated distance drops below the tolerance (useful for dwithin calculations). Return a negative distance for incalculable cases.

Definition at line 2058 of file lwgeodetic.c.

{
        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 ( FLAGS_GET_GEODETIC(lwgeom1->flags) && lwgeom1->bbox )
                gbox1 = *(lwgeom1->bbox);
        else
                lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
 
        /* Make sure we have boxes */
        if ( FLAGS_GET_GEODETIC(lwgeom2->flags) && 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;
                uint32_t 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;
                uint32_t 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)
        {
                const POINT2D* p;
                LWPOLY* lwpoly1 = (LWPOLY*)lwgeom1;
                LWPOLY* lwpoly2 = (LWPOLY*)lwgeom2;
                double distance = FLT_MAX;
                uint32_t 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) )
        {
                uint32_t 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) )
        {
                uint32_t 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;
 
}

References LWGEOM::bbox, distance(), LWGEOM::flags, FLAGS_GET_GEODETIC, gbox_init(), gbox_overlaps(), LWCOLLECTION::geoms, getPoint2d_cp(), LINETYPE, LW_FALSE, LW_TRUE, LWDEBUG, LWDEBUGF, lwerror(), lwgeom_calculate_gbox_geodetic(), lwgeom_distance_spheroid(), 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_knn(), geography_distance_uncached(), geography_dwithin(), geography_dwithin_uncached(), geometry_distance_spheroid(), lwgeom_distance_spheroid(), test_lwgeom_distance_sphere(), test_tree_circ_distance(), and test_tree_circ_distance_threshold().

Here is the call graph for this function:

Here is the caller graph for this function: