PostGIS  2.2.7dev-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 2081 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: