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 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 2092 of file lwgeodetic.c.

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: