PostGIS  2.5.0dev-r@@SVN_REVISION@@
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 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: