PostGIS  2.5.0dev-r@@SVN_REVISION@@
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 calulated distance drops below the tolerance (useful for dwithin calculations). Return a negative distance for incalculable cases.

Definition at line 2096 of file lwgeodetic.c.

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_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(), test_lwgeom_distance_sphere(), test_tree_circ_distance(), and test_tree_circ_distance_threshold().

2097 {
2098  uint8_t type1, type2;
2099  int check_intersection = LW_FALSE;
2100  GBOX gbox1, gbox2;
2101 
2102  gbox_init(&gbox1);
2103  gbox_init(&gbox2);
2104 
2105  assert(lwgeom1);
2106  assert(lwgeom2);
2107 
2108  LWDEBUGF(4, "entered function, tolerance %.8g", tolerance);
2109 
2110  /* What's the distance to an empty geometry? We don't know.
2111  Return a negative number so the caller can catch this case. */
2112  if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
2113  {
2114  return -1.0;
2115  }
2116 
2117  type1 = lwgeom1->type;
2118  type2 = lwgeom2->type;
2119 
2120  /* Make sure we have boxes */
2121  if ( lwgeom1->bbox )
2122  gbox1 = *(lwgeom1->bbox);
2123  else
2124  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
2125 
2126  /* Make sure we have boxes */
2127  if ( lwgeom2->bbox )
2128  gbox2 = *(lwgeom2->bbox);
2129  else
2130  lwgeom_calculate_gbox_geodetic(lwgeom2, &gbox2);
2131 
2132  /* If the boxes aren't disjoint, we have to check for edge intersections */
2133  if ( gbox_overlaps(&gbox1, &gbox2) )
2134  check_intersection = LW_TRUE;
2135 
2136  /* Point/line combinations can all be handled with simple point array iterations */
2137  if ( ( type1 == POINTTYPE || type1 == LINETYPE ) &&
2138  ( type2 == POINTTYPE || type2 == LINETYPE ) )
2139  {
2140  POINTARRAY *pa1, *pa2;
2141 
2142  if ( type1 == POINTTYPE )
2143  pa1 = ((LWPOINT*)lwgeom1)->point;
2144  else
2145  pa1 = ((LWLINE*)lwgeom1)->points;
2146 
2147  if ( type2 == POINTTYPE )
2148  pa2 = ((LWPOINT*)lwgeom2)->point;
2149  else
2150  pa2 = ((LWLINE*)lwgeom2)->points;
2151 
2152  return ptarray_distance_spheroid(pa1, pa2, spheroid, tolerance, check_intersection);
2153  }
2154 
2155  /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */
2156  if ( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) ||
2157  ( type2 == POLYGONTYPE && type1 == POINTTYPE ) )
2158  {
2159  const POINT2D *p;
2160  LWPOLY *lwpoly;
2161  LWPOINT *lwpt;
2162  double distance = FLT_MAX;
2163  uint32_t i;
2164 
2165  if ( type1 == POINTTYPE )
2166  {
2167  lwpt = (LWPOINT*)lwgeom1;
2168  lwpoly = (LWPOLY*)lwgeom2;
2169  }
2170  else
2171  {
2172  lwpt = (LWPOINT*)lwgeom2;
2173  lwpoly = (LWPOLY*)lwgeom1;
2174  }
2175  p = getPoint2d_cp(lwpt->point, 0);
2176 
2177  /* Point in polygon implies zero distance */
2178  if ( lwpoly_covers_point2d(lwpoly, p) )
2179  {
2180  return 0.0;
2181  }
2182 
2183  /* Not inside, so what's the actual distance? */
2184  for ( i = 0; i < lwpoly->nrings; i++ )
2185  {
2186  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwpt->point, spheroid, tolerance, check_intersection);
2187  if ( ring_distance < distance )
2188  distance = ring_distance;
2189  if ( distance < tolerance )
2190  return distance;
2191  }
2192  return distance;
2193  }
2194 
2195  /* Line/polygon case, if start point-in-poly, return zero, else return distance. */
2196  if ( ( type1 == POLYGONTYPE && type2 == LINETYPE ) ||
2197  ( type2 == POLYGONTYPE && type1 == LINETYPE ) )
2198  {
2199  const POINT2D *p;
2200  LWPOLY *lwpoly;
2201  LWLINE *lwline;
2202  double distance = FLT_MAX;
2203  uint32_t i;
2204 
2205  if ( type1 == LINETYPE )
2206  {
2207  lwline = (LWLINE*)lwgeom1;
2208  lwpoly = (LWPOLY*)lwgeom2;
2209  }
2210  else
2211  {
2212  lwline = (LWLINE*)lwgeom2;
2213  lwpoly = (LWPOLY*)lwgeom1;
2214  }
2215  p = getPoint2d_cp(lwline->points, 0);
2216 
2217  LWDEBUG(4, "checking if a point of line is in polygon");
2218 
2219  /* Point in polygon implies zero distance */
2220  if ( lwpoly_covers_point2d(lwpoly, p) )
2221  return 0.0;
2222 
2223  LWDEBUG(4, "checking ring distances");
2224 
2225  /* Not contained, so what's the actual distance? */
2226  for ( i = 0; i < lwpoly->nrings; i++ )
2227  {
2228  double ring_distance = ptarray_distance_spheroid(lwpoly->rings[i], lwline->points, spheroid, tolerance, check_intersection);
2229  LWDEBUGF(4, "ring[%d] ring_distance = %.8g", i, ring_distance);
2230  if ( ring_distance < distance )
2231  distance = ring_distance;
2232  if ( distance < tolerance )
2233  return distance;
2234  }
2235  LWDEBUGF(4, "all rings checked, returning distance = %.8g", distance);
2236  return distance;
2237 
2238  }
2239 
2240  /* Polygon/polygon case, if start point-in-poly, return zero, else
2241  * return distance. */
2242  if (type1 == POLYGONTYPE && type2 == POLYGONTYPE)
2243  {
2244  const POINT2D* p;
2245  LWPOLY* lwpoly1 = (LWPOLY*)lwgeom1;
2246  LWPOLY* lwpoly2 = (LWPOLY*)lwgeom2;
2247  double distance = FLT_MAX;
2248  uint32_t i, j;
2249 
2250  /* Point of 2 in polygon 1 implies zero distance */
2251  p = getPoint2d_cp(lwpoly1->rings[0], 0);
2252  if (lwpoly_covers_point2d(lwpoly2, p)) return 0.0;
2253 
2254  /* Point of 1 in polygon 2 implies zero distance */
2255  p = getPoint2d_cp(lwpoly2->rings[0], 0);
2256  if (lwpoly_covers_point2d(lwpoly1, p)) return 0.0;
2257 
2258  /* Not contained, so what's the actual distance? */
2259  for (i = 0; i < lwpoly1->nrings; i++)
2260  {
2261  for (j = 0; j < lwpoly2->nrings; j++)
2262  {
2263  double ring_distance =
2265  lwpoly1->rings[i],
2266  lwpoly2->rings[j],
2267  spheroid,
2268  tolerance,
2269  check_intersection);
2270  if (ring_distance < distance)
2271  distance = ring_distance;
2272  if (distance < tolerance) return distance;
2273  }
2274  }
2275  return distance;
2276  }
2277 
2278  /* Recurse into collections */
2279  if ( lwtype_is_collection(type1) )
2280  {
2281  uint32_t i;
2282  double distance = FLT_MAX;
2283  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1;
2284 
2285  for ( i = 0; i < col->ngeoms; i++ )
2286  {
2287  double geom_distance = lwgeom_distance_spheroid(
2288  col->geoms[i], lwgeom2, spheroid, tolerance);
2289  if ( geom_distance < distance )
2290  distance = geom_distance;
2291  if ( distance < tolerance )
2292  return distance;
2293  }
2294  return distance;
2295  }
2296 
2297  /* Recurse into collections */
2298  if ( lwtype_is_collection(type2) )
2299  {
2300  uint32_t i;
2301  double distance = FLT_MAX;
2302  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2;
2303 
2304  for ( i = 0; i < col->ngeoms; i++ )
2305  {
2306  double geom_distance = lwgeom_distance_spheroid(lwgeom1, col->geoms[i], spheroid, tolerance);
2307  if ( geom_distance < distance )
2308  distance = geom_distance;
2309  if ( distance < tolerance )
2310  return distance;
2311  }
2312  return distance;
2313  }
2314 
2315 
2316  lwerror("arguments include unsupported geometry type (%s, %s)", lwtype_name(type1), lwtype_name(type1));
2317  return -1.0;
2318 
2319 }
#define LINETYPE
Definition: liblwgeom.h:85
GBOX * bbox
Definition: liblwgeom.h:397
static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
Definition: lwgeodetic.c:1744
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:2096
#define POLYGONTYPE
Definition: liblwgeom.h:86
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
uint32_t ngeoms
Definition: liblwgeom.h:506
POINTARRAY * point
Definition: liblwgeom.h:410
uint32_t nrings
Definition: liblwgeom.h:454
unsigned int uint32_t
Definition: uthash.h:78
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:47
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWGEOM ** geoms
Definition: liblwgeom.h:508
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2927
POINTARRAY ** rings
Definition: liblwgeom.h:456
int lwtype_is_collection(uint8_t type)
Determine whether a type number is a collection or not.
Definition: lwgeom.c:1086
Datum distance(PG_FUNCTION_ARGS)
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
int gbox_overlaps(const GBOX *g1, const GBOX *g2)
Return LW_TRUE if the GBOX overlaps, LW_FALSE otherwise.
Definition: g_box.c:290
uint8_t type
Definition: liblwgeom.h:395
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:2427
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:1386
#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:421
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364

Here is the call graph for this function:

Here is the caller graph for this function: