PostGIS  2.1.10dev-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 2078 of file lwgeodetic.c.

References LWGEOM::bbox, distance(), gbox_init(), gbox_overlaps(), LWCOLLECTION::geoms, getPoint2d_p(), LINETYPE, LW_FALSE, LW_TRUE, LWDEBUG, LWDEBUGF, lwerror(), lwgeom_calculate_gbox_geodetic(), lwgeom_is_empty(), lwpoly_covers_point2d(), lwtype_is_collection(), lwtype_name(), MAXFLOAT, LWCOLLECTION::ngeoms, LWPOLY::nrings, LWPOINT::point, LWLINE::points, POINTTYPE, POLYGONTYPE, ptarray_distance_spheroid(), LWPOLY::rings, and LWGEOM::type.

Referenced by geography_distance(), geography_distance_uncached(), geography_dwithin(), geography_dwithin_uncached(), geometry_distance_spheroid(), test_lwgeom_distance_sphere(), and test_tree_circ_distance().

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

Here is the call graph for this function:

Here is the caller graph for this function: