PostGIS  2.4.9dev-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 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.

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