PostGIS  3.1.6dev-r@@SVN_REVISION@@

◆ 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 2190 of file lwgeodetic.c.

2191 {
2192  uint8_t type1, type2;
2193  int check_intersection = LW_FALSE;
2194  GBOX gbox1, gbox2;
2195 
2196  gbox_init(&gbox1);
2197  gbox_init(&gbox2);
2198 
2199  assert(lwgeom1);
2200  assert(lwgeom2);
2201 
2202  LWDEBUGF(4, "entered function, tolerance %.8g", tolerance);
2203 
2204  /* What's the distance to an empty geometry? We don't know.
2205  Return a negative number so the caller can catch this case. */
2206  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
2207  {
2208  return -1.0;
2209  }
2210 
2211  type1 = lwgeom1->type;
2212  type2 = lwgeom2->type;
2213 
2214  /* Make sure we have boxes */
2215  if ( lwgeom1->bbox )
2216  gbox1 = *(lwgeom1->bbox);
2217  else
2218  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
2219 
2220  /* Make sure we have boxes */
2221  if ( lwgeom2->bbox )
2222  gbox2 = *(lwgeom2->bbox);
2223  else
2224  lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);
2225 
2226  /* If the boxes aren't disjoint, we have to check for edge intersections */
2227  if ( gbox_overlaps(&gbox1, &gbox2) )
2228  check_intersection = LW_TRUE;
2229 
2230  /* Point/line combinations can all be handled with simple point array iterations */
2231  if ( ( type1 == POINTTYPE || type1 == LINETYPE ) &&
2232  ( type2 == POINTTYPE || type2 == LINETYPE ) )
2233  {
2234  POINTARRAY *pa1, *pa2;
2235 
2236  if ( type1 == POINTTYPE )
2237  pa1 = ((LWPOINT*)lwgeom1)->point;
2238  else
2239  pa1 = ((LWLINE*)lwgeom1)->points;
2240 
2241  if ( type2 == POINTTYPE )
2242  pa2 = ((LWPOINT*)lwgeom2)->point;
2243  else
2244  pa2 = ((LWLINE*)lwgeom2)->points;
2245 
2246  return ptarray_distance_spheroid(pa1, pa2, spheroid, tolerance, check_intersection);
2247  }
2248 
2249  /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */
2250  if ( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) ||
2251  ( type2 == POLYGONTYPE && type1 == POINTTYPE ) )
2252  {
2253  const POINT2D *p;
2254  LWPOLY *lwpoly;
2255  LWPOINT *lwpt;
2256  double distance = FLT_MAX;
2257  uint32_t i;
2258 
2259  if ( type1 == POINTTYPE )
2260  {
2261  lwpt = (LWPOINT*)lwgeom1;
2262  lwpoly = (LWPOLY*)lwgeom2;
2263  }
2264  else
2265  {
2266  lwpt = (LWPOINT*)lwgeom2;
2267  lwpoly = (LWPOLY*)lwgeom1;
2268  }
2269  p = getPoint2d_cp(lwpt->point, 0);
2270 
2271  /* Point in polygon implies zero distance */
2272  if ( lwpoly_covers_point2d(lwpoly, p) )
2273  {
2274  return 0.0;
2275  }
2276 
2277  /* Not inside, so what's the actual distance? */
2278  for ( i = 0; i < lwpoly->nrings; i++ )
2279  {
2280  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, spheroid, tolerance, check_intersection);
2281  if ( ring_distance < distance )
2282  distance = ring_distance;
2283  if ( distance < tolerance )
2284  return distance;
2285  }
2286  return distance;
2287  }
2288 
2289  /* Line/polygon case, if start point-in-poly, return zero, else return distance. */
2290  if ( ( type1 == POLYGONTYPE && type2 == LINETYPE ) ||
2291  ( type2 == POLYGONTYPE && type1 == LINETYPE ) )
2292  {
2293  const POINT2D *p;
2294  LWPOLY *lwpoly;
2295  LWLINE *lwline;
2296  double distance = FLT_MAX;
2297  uint32_t i;
2298 
2299  if ( type1 == LINETYPE )
2300  {
2301  lwline = (LWLINE*)lwgeom1;
2302  lwpoly = (LWPOLY*)lwgeom2;
2303  }
2304  else
2305  {
2306  lwline = (LWLINE*)lwgeom2;
2307  lwpoly = (LWPOLY*)lwgeom1;
2308  }
2309  p = getPoint2d_cp(lwline->points, 0);
2310 
2311  LWDEBUG(4, "checking if a point of line is in polygon");
2312 
2313  /* Point in polygon implies zero distance */
2314  if ( lwpoly_covers_point2d(lwpoly, p) )
2315  return 0.0;
2316 
2317  LWDEBUG(4, "checking ring distances");
2318 
2319  /* Not contained, so what's the actual distance? */
2320  for ( i = 0; i < lwpoly->nrings; i++ )
2321  {
2322  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, spheroid, tolerance, check_intersection);
2323  LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance);
2324  if ( ring_distance < distance )
2325  distance = ring_distance;
2326  if ( distance < tolerance )
2327  return distance;
2328  }
2329  LWDEBUGF(4, "all rings checked, returning distance = %.8g", distance);
2330  return distance;
2331 
2332  }
2333 
2334  /* Polygon/polygon case, if start point-in-poly, return zero, else
2335  * return distance. */
2336  if (type1 == POLYGONTYPE && type2 == POLYGONTYPE)
2337  {
2338  const POINT2D* p;
2339  LWPOLY* lwpoly1 = (LWPOLY*)lwgeom1;
2340  LWPOLY* lwpoly2 = (LWPOLY*)lwgeom2;
2341  double distance = FLT_MAX;
2342  uint32_t i, j;
2343 
2344  /* Point of 2 in polygon 1 implies zero distance */
2345  p = getPoint2d_cp(lwpoly1->rings[0], 0);
2346  if (lwpoly_covers_point2d(lwpoly2, p)) return 0.0;
2347 
2348  /* Point of 1 in polygon 2 implies zero distance */
2349  p = getPoint2d_cp(lwpoly2->rings[0], 0);
2350  if (lwpoly_covers_point2d(lwpoly1, p)) return 0.0;
2351 
2352  /* Not contained, so what's the actual distance? */
2353  for (i = 0; i < lwpoly1->nrings; i++)
2354  {
2355  for (j = 0; j < lwpoly2->nrings; j++)
2356  {
2357  double ring_distance =
2359  lwpoly1->rings[i],
2360  lwpoly2->rings[j],
2361  spheroid,
2362  tolerance,
2363  check_intersection);
2364  if (ring_distance < distance)
2365  distance = ring_distance;
2366  if (distance < tolerance) return distance;
2367  }
2368  }
2369  return distance;
2370  }
2371 
2372  /* Recurse into collections */
2373  if ( lwtype_is_collection(type1) )
2374  {
2375  uint32_t i;
2376  double distance = FLT_MAX;
2377  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
2378 
2379  for ( i = 0; i < col->ngeoms; i++ )
2380  {
2381  double geom_distance = lwgeom_distance_spheroid(
2382  col->geoms[i], lwgeom2, spheroid, tolerance);
2383  if ( geom_distance < distance )
2384  distance = geom_distance;
2385  if ( distance < tolerance )
2386  return distance;
2387  }
2388  return distance;
2389  }
2390 
2391  /* Recurse into collections */
2392  if ( lwtype_is_collection(type2) )
2393  {
2394  uint32_t i;
2395  double distance = FLT_MAX;
2396  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
2397 
2398  for ( i = 0; i < col->ngeoms; i++ )
2399  {
2400  double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], spheroid, tolerance);
2401  if ( geom_distance < distance )
2402  distance = geom_distance;
2403  if ( distance < tolerance )
2404  return distance;
2405  }
2406  return distance;
2407  }
2408 
2409 
2410  lwerror("arguments include unsupported geometry type (%s, %s)", lwtype_name(type1), lwtype_name(type1));
2411  return -1.0;
2412 
2413 }
int gbox_overlaps(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the GBOX overlaps, LW_FALSE otherwise.
Definition: gbox.c:283
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: gbox.c:40
#define LW_FALSE
Definition: liblwgeom.h:108
#define LINETYPE
Definition: liblwgeom.h:117
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition: lwgeom.c:1088
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:116
#define POLYGONTYPE
Definition: liblwgeom.h:118
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:216
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:107
int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test)
Given a polygon (lon/lat decimal degrees) and point (lon/lat decimal degrees) and a guaranteed outsid...
Definition: lwgeodetic.c:2521
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.
Definition: lwgeodetic.c:2190
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:3018
static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
Definition: lwgeodetic.c:1835
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from.
Definition: lwinline.h:101
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwinline.h:203
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
uint32_t ngeoms
Definition: liblwgeom.h:594
LWGEOM ** geoms
Definition: liblwgeom.h:589
uint8_t type
Definition: liblwgeom.h:476
GBOX * bbox
Definition: liblwgeom.h:472
POINTARRAY * points
Definition: liblwgeom.h:497
POINTARRAY * point
Definition: liblwgeom.h:485
POINTARRAY ** rings
Definition: liblwgeom.h:533
uint32_t nrings
Definition: liblwgeom.h:538

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: