PostGIS  2.5.1dev-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 2100 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().

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