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 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 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: