PostGIS  2.1.10dev-r@@SVN_REVISION@@
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 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: