PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwspheroid.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  *
5  * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
6  * Copyright (C) 2009 David Skea <David.Skea@gov.bc.ca>
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "liblwgeom_internal.h"
14 #include "lwgeodetic.h"
15 #include "lwgeom_log.h"
16 
20 void spheroid_init(SPHEROID *s, double a, double b)
21 {
22  s->a = a;
23  s->b = b;
24  s->f = (a - b) / a;
25  s->e_sq = (a*a - b*b)/(a*a);
26  s->radius = (2.0 * a + b ) / 3.0;
27 }
28 
29 static double spheroid_mu2(double alpha, const SPHEROID *s)
30 {
31  double b2 = POW2(s->b);
32  return POW2(cos(alpha)) * (POW2(s->a) - b2) / b2;
33 }
34 
35 static double spheroid_big_a(double u2)
36 {
37  return 1.0 + (u2 / 16384.0) * (4096.0 + u2 * (-768.0 + u2 * (320.0 - 175.0 * u2)));
38 }
39 
40 static double spheroid_big_b(double u2)
41 {
42  return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
43 }
44 
59 double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
60 {
61  double lambda = (b->lon - a->lon);
62  double f = spheroid->f;
63  double omf = 1 - spheroid->f;
64  double u1, u2;
65  double cos_u1, cos_u2;
66  double sin_u1, sin_u2;
67  double big_a, big_b, delta_sigma;
68  double alpha, sin_alpha, cos_alphasq, c;
69  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
70  double cos_lambda, sin_lambda;
71  double distance;
72  int i = 0;
73 
74  /* Same point => zero distance */
75  if ( geographic_point_equals(a, b) )
76  {
77  return 0.0;
78  }
79 
80  u1 = atan(omf * tan(a->lat));
81  cos_u1 = cos(u1);
82  sin_u1 = sin(u1);
83  u2 = atan(omf * tan(b->lat));
84  cos_u2 = cos(u2);
85  sin_u2 = sin(u2);
86 
87  omega = lambda;
88  do
89  {
90  cos_lambda = cos(lambda);
91  sin_lambda = sin(lambda);
92  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
93  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
94  sin_sigma = sqrt(sqrsin_sigma);
95  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
96  sigma = atan2(sin_sigma, cos_sigma);
97  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
98 
99  /* Numerical stability issue, ensure asin is not NaN */
100  if ( sin_alpha > 1.0 )
101  alpha = M_PI_2;
102  else if ( sin_alpha < -1.0 )
103  alpha = -1.0 * M_PI_2;
104  else
105  alpha = asin(sin_alpha);
106 
107  cos_alphasq = POW2(cos(alpha));
108  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
109 
110  /* Numerical stability issue, cos2 is in range */
111  if ( cos2_sigma_m > 1.0 )
112  cos2_sigma_m = 1.0;
113  if ( cos2_sigma_m < -1.0 )
114  cos2_sigma_m = -1.0;
115 
116  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
117  last_lambda = lambda;
118  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
119  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
120  i++;
121  }
122  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
123 
124  u2 = spheroid_mu2(alpha, spheroid);
125  big_a = spheroid_big_a(u2);
126  big_b = spheroid_big_b(u2);
127  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
128  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
129 
130  distance = spheroid->b * big_a * (sigma - delta_sigma);
131 
132  /* Algorithm failure, distance == NaN, fallback to sphere */
133  if ( distance != distance )
134  {
135  lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a->lat, a->lon, b->lat, b->lon, spheroid->a, spheroid->b);
136  return spheroid->radius * sphere_distance(a, b);
137  }
138 
139  return distance;
140 }
141 
155 double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
156 {
157  int i = 0;
158  double lambda = s->lon - r->lon;
159  double omf = 1 - spheroid->f;
160  double u1 = atan(omf * tan(r->lat));
161  double cos_u1 = cos(u1);
162  double sin_u1 = sin(u1);
163  double u2 = atan(omf * tan(s->lat));
164  double cos_u2 = cos(u2);
165  double sin_u2 = sin(u2);
166 
167  double omega = lambda;
168  double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
169  double sin_alpha, cos_alphasq, C, alphaFD;
170  do
171  {
172  sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
173  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
174  sin_sigma = sqrt(sqr_sin_sigma);
175  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
176  sigma = atan2(sin_sigma, cos_sigma);
177  sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
178 
179  /* Numerical stability issue, ensure asin is not NaN */
180  if ( sin_alpha > 1.0 )
181  alpha = M_PI_2;
182  else if ( sin_alpha < -1.0 )
183  alpha = -1.0 * M_PI_2;
184  else
185  alpha = asin(sin_alpha);
186 
187  cos_alphasq = POW2(cos(alpha));
188  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
189 
190  /* Numerical stability issue, cos2 is in range */
191  if ( cos2_sigma_m > 1.0 )
192  cos2_sigma_m = 1.0;
193  if ( cos2_sigma_m < -1.0 )
194  cos2_sigma_m = -1.0;
195 
196  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
197  last_lambda = lambda;
198  lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
199  (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
200  i++;
201  }
202  while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
203 
204  alphaFD = atan2((cos_u2 * sin(lambda)),
205  (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
206  if (alphaFD < 0.0)
207  {
208  alphaFD = alphaFD + 2.0 * M_PI;
209  }
210  if (alphaFD > 2.0 * M_PI)
211  {
212  alphaFD = alphaFD - 2.0 * M_PI;
213  }
214  return alphaFD;
215 }
216 
217 
232 int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
233 {
234  double omf = 1 - spheroid->f;
235  double tan_u1 = omf * tan(r->lat);
236  double u1 = atan(tan_u1);
237  double sigma, last_sigma, delta_sigma, two_sigma_m;
238  double sigma1, sin_alpha, alpha, cos_alphasq;
239  double u2, A, B;
240  double lat2, lambda, lambda2, C, omega;
241  int i = 0;
242 
243  if (azimuth < 0.0)
244  {
245  azimuth = azimuth + M_PI * 2.0;
246  }
247  if (azimuth > (PI * 2.0))
248  {
249  azimuth = azimuth - M_PI * 2.0;
250  }
251 
252  sigma1 = atan2(tan_u1, cos(azimuth));
253  sin_alpha = cos(u1) * sin(azimuth);
254  alpha = asin(sin_alpha);
255  cos_alphasq = 1.0 - POW2(sin_alpha);
256 
257  u2 = spheroid_mu2(alpha, spheroid);
258  A = spheroid_big_a(u2);
259  B = spheroid_big_b(u2);
260 
261  sigma = (distance / (spheroid->b * A));
262  do
263  {
264  two_sigma_m = 2.0 * sigma1 + sigma;
265  delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m))))));
266  last_sigma = sigma;
267  sigma = (distance / (spheroid->b * A)) + delta_sigma;
268  i++;
269  }
270  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
271 
272  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
273  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
274  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
275  cos(azimuth)))));
276  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
277  sin(u1) * sin(sigma) * cos(azimuth)));
278  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
279  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
280  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
281  lambda2 = r->lon + omega;
282  g->lat = lat2;
283  g->lon = lambda2;
284  return LW_SUCCESS;
285 }
286 
287 
288 static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
289 {
290  return spheroid->a / (sqrt(1.0 - spheroid->e_sq * POW2(sin(latitude))));
291 }
292 
293 static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
294 {
295  return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid)
296  * cos(latitude)
297  * deltaLongitude;
298 }
299 
300 
310 static double spheroid_boundingbox_area(const GEOGRAPHIC_POINT *southWestCorner, const GEOGRAPHIC_POINT *northEastCorner, const SPHEROID *spheroid)
311 {
312  double z0 = (northEastCorner->lon - southWestCorner->lon) * POW2(spheroid->b) / 2.0;
313  double e = sqrt(spheroid->e_sq);
314  double sinPhi1 = sin(southWestCorner->lat);
315  double sinPhi2 = sin(northEastCorner->lat);
316  double t1p1 = sinPhi1 / (1.0 - spheroid->e_sq * sinPhi1 * sinPhi1);
317  double t1p2 = sinPhi2 / (1.0 - spheroid->e_sq * sinPhi2 * sinPhi2);
318  double oneOver2e = 1.0 / (2.0 * e);
319  double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1));
320  double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2));
321  return z0 * (t1p2 + t2p2) - z0 * (t1p1 + t2p1);
322 }
323 
328 static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
329 {
330  GEOGRAPHIC_POINT A, B, mL, nR;
331  double deltaLng, baseArea, topArea;
332  double bE, tE, ratio, sign;
333 
334  A = *a;
335  B = *b;
336 
337  mL.lat = latitude_min;
338  mL.lon = FP_MIN(A.lon, B.lon);
339  nR.lat = FP_MIN(A.lat, B.lat);
340  nR.lon = FP_MAX(A.lon, B.lon);
341  LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
342  LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
343  baseArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
344  LWDEBUGF(4, "baseArea %.12g", baseArea);
345 
346  mL.lat = FP_MIN(A.lat, B.lat);
347  mL.lon = FP_MIN(A.lon, B.lon);
348  nR.lat = FP_MAX(A.lat, B.lat);
349  nR.lon = FP_MAX(A.lon, B.lon);
350  LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
351  LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
352  topArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
353  LWDEBUGF(4, "topArea %.12g", topArea);
354 
355  deltaLng = B.lon - A.lon;
356  LWDEBUGF(4, "deltaLng %.12g", deltaLng);
357  bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid);
358  tE = spheroid_parallel_arc_length(B.lat, deltaLng, spheroid);
359  LWDEBUGF(4, "bE %.12g", bE);
360  LWDEBUGF(4, "tE %.12g", tE);
361 
362  ratio = (bE + tE)/tE;
363  sign = signum(B.lon - A.lon);
364  return (baseArea + topArea / ratio) * sign;
365 }
366 
367 static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
368 {
369  GEOGRAPHIC_POINT a, b;
370  POINT2D p;
371  int i;
372  double area = 0.0;
373  GBOX gbox2d;
374  int in_south = LW_FALSE;
375  double delta_lon_tolerance;
376  double latitude_min;
377 
378  gbox2d.flags = gflags(0, 0, 0);
379 
380  /* Return zero on non-sensical inputs */
381  if ( ! pa || pa->npoints < 4 )
382  return 0.0;
383 
384  /* Get the raw min/max values for the latitudes */
386 
387  if ( signum(gbox2d.ymin) != signum(gbox2d.ymax) )
388  lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
389 
390  /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
391  if ( gbox2d.ymax < 0.0 )
392  in_south = LW_TRUE;
393 
394  LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
395 
396  /* Tolerance for strip area calculation */
397  if ( in_south )
398  {
399  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
400  latitude_min = deg2rad(fabs(gbox2d.ymax));
401  }
402  else
403  {
404  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
405  latitude_min = deg2rad(gbox2d.ymin);
406  }
407 
408  /* Initialize first point */
409  getPoint2d_p(pa, 0, &p);
410  geographic_point_init(p.x, p.y, &a);
411 
412  for ( i = 1; i < pa->npoints; i++ )
413  {
414  GEOGRAPHIC_POINT a1, b1;
415  double strip_area = 0.0;
416  double delta_lon = 0.0;
417  LWDEBUGF(4, "edge #%d", i);
418 
419  getPoint2d_p(pa, i, &p);
420  geographic_point_init(p.x, p.y, &b);
421 
422  a1 = a;
423  b1 = b;
424 
425  /* Flip into north if in south */
426  if ( in_south )
427  {
428  a1.lat = -1.0 * a1.lat;
429  b1.lat = -1.0 * b1.lat;
430  }
431 
432  LWDEBUGF(4, "in_south %d", in_south);
433 
434  LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
435 
436  if ( crosses_dateline(&a, &b) )
437  {
438  double shift;
439 
440  if ( a1.lon > 0.0 )
441  shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
442  else
443  shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
444 
445  LWDEBUGF(4, "shift: %.8g", shift);
446  LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
447  point_shift(&a1, shift);
448  point_shift(&b1, shift);
449  LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
450 
451  }
452 
453 
454  delta_lon = fabs(b1.lon - a1.lon);
455 
456  LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
457  LWDEBUGF(4, "delta_lon %.18g", delta_lon);
458  LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
459 
460  if ( delta_lon > 0.0 )
461  {
462  if ( delta_lon < delta_lon_tolerance )
463  {
464  strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
465  LWDEBUGF(4, "strip_area %.12g", strip_area);
466  area += strip_area;
467  }
468  else
469  {
470  GEOGRAPHIC_POINT p, q;
471  double step = floor(delta_lon / delta_lon_tolerance);
472  double distance = spheroid_distance(&a1, &b1, spheroid);
473  double pDistance = 0.0;
474  int j = 0;
475  LWDEBUGF(4, "step %.18g", step);
476  LWDEBUGF(4, "distance %.18g", distance);
477  step = distance / step;
478  LWDEBUGF(4, "step %.18g", step);
479  p = a1;
480  while (pDistance < (distance - step * 1.01))
481  {
482  double azimuth = spheroid_direction(&p, &b1, spheroid);
483  j++;
484  LWDEBUGF(4, " iteration %d", j);
485  LWDEBUGF(4, " azimuth %.12g", azimuth);
486  pDistance = pDistance + step;
487  LWDEBUGF(4, " pDistance %.12g", pDistance);
488  spheroid_project(&p, spheroid, step, azimuth, &q);
489  strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
490  LWDEBUGF(4, " strip_area %.12g", strip_area);
491  area += strip_area;
492  LWDEBUGF(4, " area %.12g", area);
493  p.lat = q.lat;
494  p.lon = q.lon;
495  }
496  strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
497  area += strip_area;
498  }
499  }
500 
501  /* B gets incremented in the next loop, so we save the value here */
502  a = b;
503  }
504  return fabs(area);
505 }
513 double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
514 {
515  int type;
516 
517  assert(lwgeom);
518 
519  /* No area in nothing */
520  if ( lwgeom_is_empty(lwgeom) )
521  return 0.0;
522 
523  /* Read the geometry type number */
524  type = lwgeom->type;
525 
526  /* Anything but polygons and collections returns zero */
527  if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
528  return 0.0;
529 
530  /* Actually calculate area */
531  if ( type == POLYGONTYPE )
532  {
533  LWPOLY *poly = (LWPOLY*)lwgeom;
534  int i;
535  double area = 0.0;
536 
537  /* Just in case there's no rings */
538  if ( poly->nrings < 1 )
539  return 0.0;
540 
541  /* First, the area of the outer ring */
542  area += ptarray_area_spheroid(poly->rings[0], spheroid);
543 
544  /* Subtract areas of inner rings */
545  for ( i = 1; i < poly->nrings; i++ )
546  {
547  area -= ptarray_area_spheroid(poly->rings[i], spheroid);
548  }
549  return area;
550  }
551 
552  /* Recurse into sub-geometries to get area */
553  if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
554  {
555  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
556  int i;
557  double area = 0.0;
558 
559  for ( i = 0; i < col->ngeoms; i++ )
560  {
561  area += lwgeom_area_spheroid(col->geoms[i], spheroid);
562  }
563  return area;
564  }
565 
566  /* Shouldn't get here. */
567  return 0.0;
568 }
569 
570 
571 
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
Calculate box (x/y) and add values to gbox.
Definition: g_box.c:471
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:40
double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
Computes the shortest distance along the surface of the spheroid between two points.
Definition: lwspheroid.c:59
double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Given two points on a unit sphere, calculate their distance apart in radians.
Definition: lwgeodetic.c:897
char * r
Definition: cu_in_wkt.c:25
int npoints
Definition: liblwgeom.h:327
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:35
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:615
int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
Given a location, an azimuth and a distance, computes the location of the projected point...
Definition: lwspheroid.c:232
#define POLYGONTYPE
Definition: liblwgeom.h:62
Datum area(PG_FUNCTION_ARGS)
static double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
Definition: lwspheroid.c:293
double b
Definition: liblwgeom.h:270
static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
This function doesn't work for edges crossing the dateline or in the southern hemisphere.
Definition: lwspheroid.c:328
#define LW_SUCCESS
Definition: liblwgeom.h:55
#define signum(a)
Ape a java function.
Definition: lwgeodetic.h:66
double radius
Definition: liblwgeom.h:274
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:29
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:137
#define FP_MIN(A, B)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:33
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double x
Definition: liblwgeom.h:284
void spheroid_init(SPHEROID *s, double a, double b)
Initialize spheroid object based on major and minor axis.
Definition: lwspheroid.c:20
double ymin
Definition: liblwgeom.h:250
double f
Definition: liblwgeom.h:271
#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
POINTARRAY ** rings
Definition: liblwgeom.h:413
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the area of an LWGEOM.
Definition: lwspheroid.c:513
int nrings
Definition: liblwgeom.h:411
double ymax
Definition: liblwgeom.h:251
double y
Definition: liblwgeom.h:284
char * s
Definition: cu_in_wkt.c:24
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
Definition: lwspheroid.c:367
double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
Computes the direction of the geodesic joining two points on the spheroid.
Definition: lwspheroid.c:155
double e_sq
Definition: liblwgeom.h:273
#define POW2(x)
Definition: lwgeodetic.h:28
Datum distance(PG_FUNCTION_ARGS)
#define PI
PI.
uint8_t flags
Definition: liblwgeom.h:247
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:131
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:65
static double spheroid_boundingbox_area(const GEOGRAPHIC_POINT *southWestCorner, const GEOGRAPHIC_POINT *northEastCorner, const SPHEROID *spheroid)
Computes the area on the spheroid of a box bounded by meridians and parallels.
Definition: lwspheroid.c:310
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:147
double a
Definition: liblwgeom.h:269
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:157
double deltaLongitude(double azimuth, double sigma, double tsm, SPHEROID *sphere)
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:60
uint8_t type
Definition: liblwgeom.h:352
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
static double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
Definition: lwspheroid.c:288
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
#define FP_MAX(A, B)