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