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

Definition at line 2092 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_distance(), 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().

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