PostGIS  2.2.7dev-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 
17 /* GeographicLib */
18 #if PROJ_GEODESIC
19 #include <geodesic.h>
20 #endif
21 
25 void spheroid_init(SPHEROID *s, double a, double b)
26 {
27  s->a = a;
28  s->b = b;
29  s->f = (a - b) / a;
30  s->e_sq = (a*a - b*b)/(a*a);
31  s->radius = (2.0 * a + b ) / 3.0;
32 }
33 
34 #if ! PROJ_GEODESIC
35 static double spheroid_mu2(double alpha, const SPHEROID *s)
36 {
37  double b2 = POW2(s->b);
38  return POW2(cos(alpha)) * (POW2(s->a) - b2) / b2;
39 }
40 
41 static double spheroid_big_a(double u2)
42 {
43  return 1.0 + (u2 / 16384.0) * (4096.0 + u2 * (-768.0 + u2 * (320.0 - 175.0 * u2)));
44 }
45 
46 static double spheroid_big_b(double u2)
47 {
48  return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
49 }
50 #endif /* ! PROJ_GEODESIC */
51 
52 
53 #if PROJ_GEODESIC
54 
65 double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
66 {
67  struct geod_geodesic gd;
68  geod_init(&gd, spheroid->a, spheroid->f);
69  double lat1 = a->lat * 180.0 / M_PI;
70  double lon1 = a->lon * 180.0 / M_PI;
71  double lat2 = b->lat * 180.0 / M_PI;
72  double lon2 = b->lon * 180.0 / M_PI;
73  double s12; /* return distance */
74  geod_inverse(&gd, lat1, lon1, lat2, lon2, &s12, 0, 0);
75  return s12;
76 }
77 
86 double spheroid_direction(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
87 {
88  struct geod_geodesic gd;
89  geod_init(&gd, spheroid->a, spheroid->f);
90  double lat1 = a->lat * 180.0 / M_PI;
91  double lon1 = a->lon * 180.0 / M_PI;
92  double lat2 = b->lat * 180.0 / M_PI;
93  double lon2 = b->lon * 180.0 / M_PI;
94  double azi1; /* return azimuth */
95  geod_inverse(&gd, lat1, lon1, lat2, lon2, 0, &azi1, 0);
96  return azi1 * M_PI / 180.0;
97 }
98 
109 int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
110 {
111  struct geod_geodesic gd;
112  geod_init(&gd, spheroid->a, spheroid->f);
113  double lat1 = r->lat * 180.0 / M_PI;
114  double lon1 = r->lon * 180.0 / M_PI;
115  double lat2, lon2; /* return projected position */
116  geod_direct(&gd, lat1, lon1, azimuth * 180.0 / M_PI, distance, &lat2, &lon2, 0);
117  g->lat = lat2 * M_PI / 180.0;
118  g->lon = lon2 * M_PI / 180.0;
119  return LW_SUCCESS;
120 }
121 
122 
123 static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
124 {
125  /* Return zero on non-sensical inputs */
126  if ( ! pa || pa->npoints < 4 )
127  return 0.0;
128 
129  struct geod_geodesic gd;
130  geod_init(&gd, spheroid->a, spheroid->f);
131  struct geod_polygon poly;
132  geod_polygon_init(&poly, 0);
133  int i;
134  double area; /* returned polygon area */
135  POINT2D p; /* long/lat units are degrees */
136 
137  /* Pass points from point array; don't close the linearring */
138  for ( i = 0; i < pa->npoints - 1; i++ )
139  {
140  getPoint2d_p(pa, i, &p);
141  geod_polygon_addpoint(&gd, &poly, p.y, p.x);
142  LWDEBUGF(4, "geod_polygon_addpoint %d: %.12g %.12g", i, p.y, p.x);
143  }
144  i = geod_polygon_compute(&gd, &poly, 0, 1, &area, 0);
145  if ( i != pa->npoints - 1 )
146  {
147  lwerror("ptarray_area_spheroid: different number of points %d vs %d",
148  i, pa->npoints - 1);
149  }
150  LWDEBUGF(4, "geod_polygon_compute area: %.12g", area);
151  return fabs(area);
152 }
153 
154 /* Above use GeographicLib */
155 #else /* ! PROJ_GEODESIC */
156 /* Below use pre-version 2.2 geodesic functions */
157 
172 double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
173 {
174  double lambda = (b->lon - a->lon);
175  double f = spheroid->f;
176  double omf = 1 - spheroid->f;
177  double u1, u2;
178  double cos_u1, cos_u2;
179  double sin_u1, sin_u2;
180  double big_a, big_b, delta_sigma;
181  double alpha, sin_alpha, cos_alphasq, c;
182  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
183  double cos_lambda, sin_lambda;
184  double distance;
185  int i = 0;
186 
187  /* Same point => zero distance */
188  if ( geographic_point_equals(a, b) )
189  {
190  return 0.0;
191  }
192 
193  u1 = atan(omf * tan(a->lat));
194  cos_u1 = cos(u1);
195  sin_u1 = sin(u1);
196  u2 = atan(omf * tan(b->lat));
197  cos_u2 = cos(u2);
198  sin_u2 = sin(u2);
199 
200  omega = lambda;
201  do
202  {
203  cos_lambda = cos(lambda);
204  sin_lambda = sin(lambda);
205  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
206  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
207  sin_sigma = sqrt(sqrsin_sigma);
208  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
209  sigma = atan2(sin_sigma, cos_sigma);
210  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
211 
212  /* Numerical stability issue, ensure asin is not NaN */
213  if ( sin_alpha > 1.0 )
214  alpha = M_PI_2;
215  else if ( sin_alpha < -1.0 )
216  alpha = -1.0 * M_PI_2;
217  else
218  alpha = asin(sin_alpha);
219 
220  cos_alphasq = POW2(cos(alpha));
221  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
222 
223  /* Numerical stability issue, cos2 is in range */
224  if ( cos2_sigma_m > 1.0 )
225  cos2_sigma_m = 1.0;
226  if ( cos2_sigma_m < -1.0 )
227  cos2_sigma_m = -1.0;
228 
229  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
230  last_lambda = lambda;
231  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
232  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
233  i++;
234  }
235  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
236 
237  u2 = spheroid_mu2(alpha, spheroid);
238  big_a = spheroid_big_a(u2);
239  big_b = spheroid_big_b(u2);
240  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
241  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
242 
243  distance = spheroid->b * big_a * (sigma - delta_sigma);
244 
245  /* Algorithm failure, distance == NaN, fallback to sphere */
246  if ( distance != distance )
247  {
248  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);
249  return spheroid->radius * sphere_distance(a, b);
250  }
251 
252  return distance;
253 }
254 
268 double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
269 {
270  int i = 0;
271  double lambda = s->lon - r->lon;
272  double omf = 1 - spheroid->f;
273  double u1 = atan(omf * tan(r->lat));
274  double cos_u1 = cos(u1);
275  double sin_u1 = sin(u1);
276  double u2 = atan(omf * tan(s->lat));
277  double cos_u2 = cos(u2);
278  double sin_u2 = sin(u2);
279 
280  double omega = lambda;
281  double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
282  double sin_alpha, cos_alphasq, C, alphaFD;
283  do
284  {
285  sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
286  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
287  sin_sigma = sqrt(sqr_sin_sigma);
288  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
289  sigma = atan2(sin_sigma, cos_sigma);
290  sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
291 
292  /* Numerical stability issue, ensure asin is not NaN */
293  if ( sin_alpha > 1.0 )
294  alpha = M_PI_2;
295  else if ( sin_alpha < -1.0 )
296  alpha = -1.0 * M_PI_2;
297  else
298  alpha = asin(sin_alpha);
299 
300  cos_alphasq = POW2(cos(alpha));
301  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
302 
303  /* Numerical stability issue, cos2 is in range */
304  if ( cos2_sigma_m > 1.0 )
305  cos2_sigma_m = 1.0;
306  if ( cos2_sigma_m < -1.0 )
307  cos2_sigma_m = -1.0;
308 
309  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
310  last_lambda = lambda;
311  lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
312  (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
313  i++;
314  }
315  while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
316 
317  alphaFD = atan2((cos_u2 * sin(lambda)),
318  (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
319  if (alphaFD < 0.0)
320  {
321  alphaFD = alphaFD + 2.0 * M_PI;
322  }
323  if (alphaFD > 2.0 * M_PI)
324  {
325  alphaFD = alphaFD - 2.0 * M_PI;
326  }
327  return alphaFD;
328 }
329 
330 
345 int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
346 {
347  double omf = 1 - spheroid->f;
348  double tan_u1 = omf * tan(r->lat);
349  double u1 = atan(tan_u1);
350  double sigma, last_sigma, delta_sigma, two_sigma_m;
351  double sigma1, sin_alpha, alpha, cos_alphasq;
352  double u2, A, B;
353  double lat2, lambda, lambda2, C, omega;
354  int i = 0;
355 
356  if (azimuth < 0.0)
357  {
358  azimuth = azimuth + M_PI * 2.0;
359  }
360  if (azimuth > (M_PI * 2.0))
361  {
362  azimuth = azimuth - M_PI * 2.0;
363  }
364 
365  sigma1 = atan2(tan_u1, cos(azimuth));
366  sin_alpha = cos(u1) * sin(azimuth);
367  alpha = asin(sin_alpha);
368  cos_alphasq = 1.0 - POW2(sin_alpha);
369 
370  u2 = spheroid_mu2(alpha, spheroid);
371  A = spheroid_big_a(u2);
372  B = spheroid_big_b(u2);
373 
374  sigma = (distance / (spheroid->b * A));
375  do
376  {
377  two_sigma_m = 2.0 * sigma1 + sigma;
378  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))))));
379  last_sigma = sigma;
380  sigma = (distance / (spheroid->b * A)) + delta_sigma;
381  i++;
382  }
383  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
384 
385  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
386  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
387  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
388  cos(azimuth)))));
389  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
390  sin(u1) * sin(sigma) * cos(azimuth)));
391  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
392  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
393  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
394  lambda2 = r->lon + omega;
395  g->lat = lat2;
396  g->lon = lambda2;
397  return LW_SUCCESS;
398 }
399 
400 
401 static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
402 {
403  return spheroid->a / (sqrt(1.0 - spheroid->e_sq * POW2(sin(latitude))));
404 }
405 
406 static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
407 {
408  return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid)
409  * cos(latitude)
410  * deltaLongitude;
411 }
412 
413 
423 static double spheroid_boundingbox_area(const GEOGRAPHIC_POINT *southWestCorner, const GEOGRAPHIC_POINT *northEastCorner, const SPHEROID *spheroid)
424 {
425  double z0 = (northEastCorner->lon - southWestCorner->lon) * POW2(spheroid->b) / 2.0;
426  double e = sqrt(spheroid->e_sq);
427  double sinPhi1 = sin(southWestCorner->lat);
428  double sinPhi2 = sin(northEastCorner->lat);
429  double t1p1 = sinPhi1 / (1.0 - spheroid->e_sq * sinPhi1 * sinPhi1);
430  double t1p2 = sinPhi2 / (1.0 - spheroid->e_sq * sinPhi2 * sinPhi2);
431  double oneOver2e = 1.0 / (2.0 * e);
432  double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1));
433  double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2));
434  return z0 * (t1p2 + t2p2) - z0 * (t1p1 + t2p1);
435 }
436 
441 static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
442 {
443  GEOGRAPHIC_POINT A, B, mL, nR;
444  double deltaLng, baseArea, topArea;
445  double bE, tE, ratio, sign;
446 
447  A = *a;
448  B = *b;
449 
450  mL.lat = latitude_min;
451  mL.lon = FP_MIN(A.lon, B.lon);
452  nR.lat = FP_MIN(A.lat, B.lat);
453  nR.lon = FP_MAX(A.lon, B.lon);
454  LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
455  LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
456  baseArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
457  LWDEBUGF(4, "baseArea %.12g", baseArea);
458 
459  mL.lat = FP_MIN(A.lat, B.lat);
460  mL.lon = FP_MIN(A.lon, B.lon);
461  nR.lat = FP_MAX(A.lat, B.lat);
462  nR.lon = FP_MAX(A.lon, B.lon);
463  LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
464  LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
465  topArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
466  LWDEBUGF(4, "topArea %.12g", topArea);
467 
468  deltaLng = B.lon - A.lon;
469  LWDEBUGF(4, "deltaLng %.12g", deltaLng);
470  bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid);
471  tE = spheroid_parallel_arc_length(B.lat, deltaLng, spheroid);
472  LWDEBUGF(4, "bE %.12g", bE);
473  LWDEBUGF(4, "tE %.12g", tE);
474 
475  ratio = (bE + tE)/tE;
476  sign = SIGNUM(B.lon - A.lon);
477  return (baseArea + topArea / ratio) * sign;
478 }
479 
480 static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
481 {
482  GEOGRAPHIC_POINT a, b;
483  POINT2D p;
484  int i;
485  double area = 0.0;
486  GBOX gbox2d;
487  int in_south = LW_FALSE;
488  double delta_lon_tolerance;
489  double latitude_min;
490 
491  gbox2d.flags = gflags(0, 0, 0);
492 
493  /* Return zero on non-sensical inputs */
494  if ( ! pa || pa->npoints < 4 )
495  return 0.0;
496 
497  /* Get the raw min/max values for the latitudes */
499 
500  if ( SIGNUM(gbox2d.ymin) != SIGNUM(gbox2d.ymax) )
501  lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
502 
503  /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
504  if ( gbox2d.ymax < 0.0 )
505  in_south = LW_TRUE;
506 
507  LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
508 
509  /* Tolerance for strip area calculation */
510  if ( in_south )
511  {
512  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
513  latitude_min = deg2rad(fabs(gbox2d.ymax));
514  }
515  else
516  {
517  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
518  latitude_min = deg2rad(gbox2d.ymin);
519  }
520 
521  /* Initialize first point */
522  getPoint2d_p(pa, 0, &p);
523  geographic_point_init(p.x, p.y, &a);
524 
525  for ( i = 1; i < pa->npoints; i++ )
526  {
527  GEOGRAPHIC_POINT a1, b1;
528  double strip_area = 0.0;
529  double delta_lon = 0.0;
530  LWDEBUGF(4, "edge #%d", i);
531 
532  getPoint2d_p(pa, i, &p);
533  geographic_point_init(p.x, p.y, &b);
534 
535  a1 = a;
536  b1 = b;
537 
538  /* Flip into north if in south */
539  if ( in_south )
540  {
541  a1.lat = -1.0 * a1.lat;
542  b1.lat = -1.0 * b1.lat;
543  }
544 
545  LWDEBUGF(4, "in_south %d", in_south);
546 
547  LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
548 
549  if ( crosses_dateline(&a, &b) )
550  {
551  double shift;
552 
553  if ( a1.lon > 0.0 )
554  shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
555  else
556  shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
557 
558  LWDEBUGF(4, "shift: %.8g", shift);
559  LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
560  point_shift(&a1, shift);
561  point_shift(&b1, shift);
562  LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
563 
564  }
565 
566 
567  delta_lon = fabs(b1.lon - a1.lon);
568 
569  LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
570  LWDEBUGF(4, "delta_lon %.18g", delta_lon);
571  LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
572 
573  if ( delta_lon > 0.0 )
574  {
575  if ( delta_lon < delta_lon_tolerance )
576  {
577  strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
578  LWDEBUGF(4, "strip_area %.12g", strip_area);
579  area += strip_area;
580  }
581  else
582  {
583  GEOGRAPHIC_POINT p, q;
584  double step = floor(delta_lon / delta_lon_tolerance);
585  double distance = spheroid_distance(&a1, &b1, spheroid);
586  double pDistance = 0.0;
587  int j = 0;
588  LWDEBUGF(4, "step %.18g", step);
589  LWDEBUGF(4, "distance %.18g", distance);
590  step = distance / step;
591  LWDEBUGF(4, "step %.18g", step);
592  p = a1;
593  while (pDistance < (distance - step * 1.01))
594  {
595  double azimuth = spheroid_direction(&p, &b1, spheroid);
596  j++;
597  LWDEBUGF(4, " iteration %d", j);
598  LWDEBUGF(4, " azimuth %.12g", azimuth);
599  pDistance = pDistance + step;
600  LWDEBUGF(4, " pDistance %.12g", pDistance);
601  spheroid_project(&p, spheroid, step, azimuth, &q);
602  strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
603  LWDEBUGF(4, " strip_area %.12g", strip_area);
604  area += strip_area;
605  LWDEBUGF(4, " area %.12g", area);
606  p.lat = q.lat;
607  p.lon = q.lon;
608  }
609  strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
610  area += strip_area;
611  }
612  }
613 
614  /* B gets incremented in the next loop, so we save the value here */
615  a = b;
616  }
617  return fabs(area);
618 }
619 #endif /* else ! PROJ_GEODESIC */
620 
628 double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
629 {
630  int type;
631 
632  assert(lwgeom);
633 
634  /* No area in nothing */
635  if ( lwgeom_is_empty(lwgeom) )
636  return 0.0;
637 
638  /* Read the geometry type number */
639  type = lwgeom->type;
640 
641  /* Anything but polygons and collections returns zero */
642  if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
643  return 0.0;
644 
645  /* Actually calculate area */
646  if ( type == POLYGONTYPE )
647  {
648  LWPOLY *poly = (LWPOLY*)lwgeom;
649  int i;
650  double area = 0.0;
651 
652  /* Just in case there's no rings */
653  if ( poly->nrings < 1 )
654  return 0.0;
655 
656  /* First, the area of the outer ring */
657  area += ptarray_area_spheroid(poly->rings[0], spheroid);
658 
659  /* Subtract areas of inner rings */
660  for ( i = 1; i < poly->nrings; i++ )
661  {
662  area -= ptarray_area_spheroid(poly->rings[i], spheroid);
663  }
664  return area;
665  }
666 
667  /* Recurse into sub-geometries to get area */
668  if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
669  {
670  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
671  int i;
672  double area = 0.0;
673 
674  for ( i = 0; i < col->ngeoms; i++ )
675  {
676  area += lwgeom_area_spheroid(col->geoms[i], spheroid);
677  }
678  return area;
679  }
680 
681  /* Shouldn't get here. */
682  return 0.0;
683 }
684 
685 
686 
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
Calculate box (x/y) and add values to gbox.
Definition: g_box.c:512
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:46
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:172
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:898
char * r
Definition: cu_in_wkt.c:24
int npoints
Definition: liblwgeom.h:355
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:41
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:616
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:345
#define POLYGONTYPE
Definition: liblwgeom.h:72
Datum area(PG_FUNCTION_ARGS)
static double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
Definition: lwspheroid.c:406
double b
Definition: liblwgeom.h:298
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:441
#define LW_SUCCESS
Definition: liblwgeom.h:65
double radius
Definition: liblwgeom.h:302
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:35
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:136
#define FP_MIN(A, B)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:32
double x
Definition: liblwgeom.h:312
void spheroid_init(SPHEROID *s, double a, double b)
Initialize spheroid object based on major and minor axis.
Definition: lwspheroid.c:25
double ymin
Definition: liblwgeom.h:278
double f
Definition: liblwgeom.h:299
#define LW_FALSE
Definition: liblwgeom.h:62
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
LWGEOM ** geoms
Definition: liblwgeom.h:493
POINTARRAY ** rings
Definition: liblwgeom.h:441
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the area of an LWGEOM.
Definition: lwspheroid.c:628
int nrings
Definition: liblwgeom.h:439
double ymax
Definition: liblwgeom.h:279
double y
Definition: liblwgeom.h:312
char * s
Definition: cu_in_wkt.c:23
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:448
static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
Definition: lwspheroid.c:480
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:268
double e_sq
Definition: liblwgeom.h:301
#define POW2(x)
Definition: lwgeodetic.h:27
Datum distance(PG_FUNCTION_ARGS)
uint8_t flags
Definition: liblwgeom.h:275
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:130
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:75
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:423
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:146
double a
Definition: liblwgeom.h:297
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:156
double deltaLongitude(double azimuth, double sigma, double tsm, SPHEROID *sphere)
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:59
uint8_t type
Definition: liblwgeom.h:380
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
static double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
Definition: lwspheroid.c:401
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:74
#define COLLECTIONTYPE
Definition: liblwgeom.h:76
#define SIGNUM(n)
Macro that returns: -1 if n < 0, 1 if n > 0, 0 if n == 0.
#define FP_MAX(A, B)