lwgeom_distance_spheroid()

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

Calculate the distance between two LWGEOMs, using the coordinates are longitude and latitude.

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 2204 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 ( 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;
                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(), 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: