PostGIS  2.5.7dev-r@@SVN_REVISION@@
lwspheroid.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
22  * Copyright (C) 2009 David Skea <David.Skea@gov.bc.ca>
23  *
24  **********************************************************************/
25 
26 
27 #include "liblwgeom_internal.h"
28 #include "lwgeodetic.h"
29 #include "lwgeom_log.h"
30 
31 /* GeographicLib */
32 #if PROJ_GEODESIC
33 #include <geodesic.h>
34 #endif
35 
39 void spheroid_init(SPHEROID *s, double a, double b)
40 {
41  s->a = a;
42  s->b = b;
43  s->f = (a - b) / a;
44  s->e_sq = (a*a - b*b)/(a*a);
45  s->radius = (2.0 * a + b ) / 3.0;
46 }
47 
48 #if ! PROJ_GEODESIC
49 static double spheroid_mu2(double alpha, const SPHEROID *s)
50 {
51  double b2 = POW2(s->b);
52  return POW2(cos(alpha)) * (POW2(s->a) - b2) / b2;
53 }
54 
55 static double spheroid_big_a(double u2)
56 {
57  return 1.0 + (u2 / 16384.0) * (4096.0 + u2 * (-768.0 + u2 * (320.0 - 175.0 * u2)));
58 }
59 
60 static double spheroid_big_b(double u2)
61 {
62  return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
63 }
64 #endif /* ! PROJ_GEODESIC */
65 
66 
67 #if PROJ_GEODESIC
68 
79 double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
80 {
81  struct geod_geodesic gd;
82 
83  /* Same point => zero distance */
84  if ( geographic_point_equals(a, b) )
85  return 0.0;
86 
87  geod_init(&gd, spheroid->a, spheroid->f);
88  double lat1 = a->lat * 180.0 / M_PI;
89  double lon1 = a->lon * 180.0 / M_PI;
90  double lat2 = b->lat * 180.0 / M_PI;
91  double lon2 = b->lon * 180.0 / M_PI;
92  double s12 = 0.0; /* return distance */
93  geod_inverse(&gd, lat1, lon1, lat2, lon2, &s12, 0, 0);
94  return s12;
95 }
96 
105 double spheroid_direction(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
106 {
107  struct geod_geodesic gd;
108  geod_init(&gd, spheroid->a, spheroid->f);
109  double lat1 = a->lat * 180.0 / M_PI;
110  double lon1 = a->lon * 180.0 / M_PI;
111  double lat2 = b->lat * 180.0 / M_PI;
112  double lon2 = b->lon * 180.0 / M_PI;
113  double azi1; /* return azimuth */
114  geod_inverse(&gd, lat1, lon1, lat2, lon2, 0, &azi1, 0);
115  return azi1 * M_PI / 180.0;
116 }
117 
128 int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
129 {
130  struct geod_geodesic gd;
131  geod_init(&gd, spheroid->a, spheroid->f);
132  double lat1 = r->lat * 180.0 / M_PI;
133  double lon1 = r->lon * 180.0 / M_PI;
134  double lat2, lon2; /* return projected position */
135  geod_direct(&gd, lat1, lon1, azimuth * 180.0 / M_PI, distance, &lat2, &lon2, 0);
136  g->lat = lat2 * M_PI / 180.0;
137  g->lon = lon2 * M_PI / 180.0;
138  return LW_SUCCESS;
139 }
140 
141 
142 static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
143 {
144  /* Return zero on non-sensical inputs */
145  if ( ! pa || pa->npoints < 4 )
146  return 0.0;
147 
148  struct geod_geodesic gd;
149  geod_init(&gd, spheroid->a, spheroid->f);
150  struct geod_polygon poly;
151  geod_polygon_init(&poly, 0);
152  uint32_t i;
153  double area; /* returned polygon area */
154  POINT2D p; /* long/lat units are degrees */
155 
156  /* Pass points from point array; don't close the linearring */
157  for ( i = 0; i < pa->npoints - 1; i++ )
158  {
159  getPoint2d_p(pa, i, &p);
160  geod_polygon_addpoint(&gd, &poly, p.y, p.x);
161  LWDEBUGF(4, "geod_polygon_addpoint %d: %.12g %.12g", i, p.y, p.x);
162  }
163  i = geod_polygon_compute(&gd, &poly, 0, 1, &area, 0);
164  if ( i != pa->npoints - 1 )
165  {
166  lwerror("ptarray_area_spheroid: different number of points %d vs %d",
167  i, pa->npoints - 1);
168  }
169  LWDEBUGF(4, "geod_polygon_compute area: %.12g", area);
170  return fabs(area);
171 }
172 
173 /* Above use GeographicLib */
174 #else /* ! PROJ_GEODESIC */
175 /* Below use pre-version 2.2 geodesic functions */
176 
191 double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
192 {
193  double lambda = (b->lon - a->lon);
194  double f = spheroid->f;
195  double omf = 1 - spheroid->f;
196  double u1, u2;
197  double cos_u1, cos_u2;
198  double sin_u1, sin_u2;
199  double big_a, big_b, delta_sigma;
200  double alpha, sin_alpha, cos_alphasq, c;
201  double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
202  double cos_lambda, sin_lambda;
203  double distance;
204  int i = 0;
205 
206  /* Same point => zero distance */
207  if ( geographic_point_equals(a, b) )
208  {
209  return 0.0;
210  }
211 
212  u1 = atan(omf * tan(a->lat));
213  cos_u1 = cos(u1);
214  sin_u1 = sin(u1);
215  u2 = atan(omf * tan(b->lat));
216  cos_u2 = cos(u2);
217  sin_u2 = sin(u2);
218 
219  omega = lambda;
220  do
221  {
222  cos_lambda = cos(lambda);
223  sin_lambda = sin(lambda);
224  sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
225  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
226  sin_sigma = sqrt(sqrsin_sigma);
227  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
228  sigma = atan2(sin_sigma, cos_sigma);
229  sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
230 
231  /* Numerical stability issue, ensure asin is not NaN */
232  if ( sin_alpha > 1.0 )
233  alpha = M_PI_2;
234  else if ( sin_alpha < -1.0 )
235  alpha = -1.0 * M_PI_2;
236  else
237  alpha = asin(sin_alpha);
238 
239  cos_alphasq = POW2(cos(alpha));
240  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
241 
242  /* Numerical stability issue, cos2 is in range */
243  if ( cos2_sigma_m > 1.0 )
244  cos2_sigma_m = 1.0;
245  if ( cos2_sigma_m < -1.0 )
246  cos2_sigma_m = -1.0;
247 
248  c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
249  last_lambda = lambda;
250  lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
251  (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
252  i++;
253  }
254  while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
255 
256  u2 = spheroid_mu2(alpha, spheroid);
257  big_a = spheroid_big_a(u2);
258  big_b = spheroid_big_b(u2);
259  delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
260  (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
261 
262  distance = spheroid->b * big_a * (sigma - delta_sigma);
263 
264  /* Algorithm failure, distance == NaN, fallback to sphere */
265  if ( distance != distance )
266  {
267  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);
268  return spheroid->radius * sphere_distance(a, b);
269  }
270 
271  return distance;
272 }
273 
287 double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
288 {
289  int i = 0;
290  double lambda = s->lon - r->lon;
291  double omf = 1 - spheroid->f;
292  double u1 = atan(omf * tan(r->lat));
293  double cos_u1 = cos(u1);
294  double sin_u1 = sin(u1);
295  double u2 = atan(omf * tan(s->lat));
296  double cos_u2 = cos(u2);
297  double sin_u2 = sin(u2);
298 
299  double omega = lambda;
300  double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
301  double sin_alpha, cos_alphasq, C, alphaFD;
302  do
303  {
304  sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
305  POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
306  sin_sigma = sqrt(sqr_sin_sigma);
307  cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
308  sigma = atan2(sin_sigma, cos_sigma);
309  sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
310 
311  /* Numerical stability issue, ensure asin is not NaN */
312  if ( sin_alpha > 1.0 )
313  alpha = M_PI_2;
314  else if ( sin_alpha < -1.0 )
315  alpha = -1.0 * M_PI_2;
316  else
317  alpha = asin(sin_alpha);
318 
319  cos_alphasq = POW2(cos(alpha));
320  cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
321 
322  /* Numerical stability issue, cos2 is in range */
323  if ( cos2_sigma_m > 1.0 )
324  cos2_sigma_m = 1.0;
325  if ( cos2_sigma_m < -1.0 )
326  cos2_sigma_m = -1.0;
327 
328  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
329  last_lambda = lambda;
330  lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
331  (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
332  i++;
333  }
334  while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
335 
336  alphaFD = atan2((cos_u2 * sin(lambda)),
337  (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
338  if (alphaFD < 0.0)
339  {
340  alphaFD = alphaFD + 2.0 * M_PI;
341  }
342  if (alphaFD > 2.0 * M_PI)
343  {
344  alphaFD = alphaFD - 2.0 * M_PI;
345  }
346  return alphaFD;
347 }
348 
349 
364 int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
365 {
366  double omf = 1 - spheroid->f;
367  double tan_u1 = omf * tan(r->lat);
368  double u1 = atan(tan_u1);
369  double sigma, last_sigma, delta_sigma, two_sigma_m;
370  double sigma1, sin_alpha, alpha, cos_alphasq;
371  double u2, A, B;
372  double lat2, lambda, lambda2, C, omega;
373  int i = 0;
374 
375  if (azimuth < 0.0)
376  {
377  azimuth = azimuth + M_PI * 2.0;
378  }
379  if (azimuth > (M_PI * 2.0))
380  {
381  azimuth = azimuth - M_PI * 2.0;
382  }
383 
384  sigma1 = atan2(tan_u1, cos(azimuth));
385  sin_alpha = cos(u1) * sin(azimuth);
386  alpha = asin(sin_alpha);
387  cos_alphasq = 1.0 - POW2(sin_alpha);
388 
389  u2 = spheroid_mu2(alpha, spheroid);
390  A = spheroid_big_a(u2);
391  B = spheroid_big_b(u2);
392 
393  sigma = (distance / (spheroid->b * A));
394  do
395  {
396  two_sigma_m = 2.0 * sigma1 + sigma;
397  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))))));
398  last_sigma = sigma;
399  sigma = (distance / (spheroid->b * A)) + delta_sigma;
400  i++;
401  }
402  while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
403 
404  lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
405  cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
406  POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
407  cos(azimuth)))));
408  lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
409  sin(u1) * sin(sigma) * cos(azimuth)));
410  C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
411  omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
412  (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
413  lambda2 = r->lon + omega;
414  g->lat = lat2;
415  g->lon = lambda2;
416  return LW_SUCCESS;
417 }
418 
419 
420 static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
421 {
422  return spheroid->a / (sqrt(1.0 - spheroid->e_sq * POW2(sin(latitude))));
423 }
424 
425 static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
426 {
427  return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid)
428  * cos(latitude)
429  * deltaLongitude;
430 }
431 
432 
442 static double spheroid_boundingbox_area(const GEOGRAPHIC_POINT *southWestCorner, const GEOGRAPHIC_POINT *northEastCorner, const SPHEROID *spheroid)
443 {
444  double z0 = (northEastCorner->lon - southWestCorner->lon) * POW2(spheroid->b) / 2.0;
445  double e = sqrt(spheroid->e_sq);
446  double sinPhi1 = sin(southWestCorner->lat);
447  double sinPhi2 = sin(northEastCorner->lat);
448  double t1p1 = sinPhi1 / (1.0 - spheroid->e_sq * sinPhi1 * sinPhi1);
449  double t1p2 = sinPhi2 / (1.0 - spheroid->e_sq * sinPhi2 * sinPhi2);
450  double oneOver2e = 1.0 / (2.0 * e);
451  double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1));
452  double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2));
453  return z0 * (t1p2 + t2p2) - z0 * (t1p1 + t2p1);
454 }
455 
460 static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
461 {
462  GEOGRAPHIC_POINT A, B, mL, nR;
463  double deltaLng, baseArea, topArea;
464  double bE, tE, ratio, sign;
465 
466  A = *a;
467  B = *b;
468 
469  mL.lat = latitude_min;
470  mL.lon = FP_MIN(A.lon, B.lon);
471  nR.lat = FP_MIN(A.lat, B.lat);
472  nR.lon = FP_MAX(A.lon, B.lon);
473  LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
474  LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
475  baseArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
476  LWDEBUGF(4, "baseArea %.12g", baseArea);
477 
478  mL.lat = FP_MIN(A.lat, B.lat);
479  mL.lon = FP_MIN(A.lon, B.lon);
480  nR.lat = FP_MAX(A.lat, B.lat);
481  nR.lon = FP_MAX(A.lon, B.lon);
482  LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
483  LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
484  topArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
485  LWDEBUGF(4, "topArea %.12g", topArea);
486 
487  deltaLng = B.lon - A.lon;
488  LWDEBUGF(4, "deltaLng %.12g", deltaLng);
489  bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid);
490  tE = spheroid_parallel_arc_length(B.lat, deltaLng, spheroid);
491  LWDEBUGF(4, "bE %.12g", bE);
492  LWDEBUGF(4, "tE %.12g", tE);
493 
494  ratio = (bE + tE)/tE;
495  sign = SIGNUM(B.lon - A.lon);
496  return (baseArea + topArea / ratio) * sign;
497 }
498 
499 static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
500 {
501  GEOGRAPHIC_POINT a, b;
502  POINT2D p;
503  int i;
504  double area = 0.0;
505  GBOX gbox2d;
506  int in_south = LW_FALSE;
507  double delta_lon_tolerance;
508  double latitude_min;
509 
510  gbox2d.flags = gflags(0, 0, 0);
511 
512  /* Return zero on non-sensical inputs */
513  if ( ! pa || pa->npoints < 4 )
514  return 0.0;
515 
516  /* Get the raw min/max values for the latitudes */
518 
519  if ( SIGNUM(gbox2d.ymin) != SIGNUM(gbox2d.ymax) )
520  lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
521 
522  /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
523  if ( gbox2d.ymax < 0.0 )
524  in_south = LW_TRUE;
525 
526  LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
527 
528  /* Tolerance for strip area calculation */
529  if ( in_south )
530  {
531  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
532  latitude_min = deg2rad(fabs(gbox2d.ymax));
533  }
534  else
535  {
536  delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
537  latitude_min = deg2rad(gbox2d.ymin);
538  }
539 
540  /* Initialize first point */
541  getPoint2d_p(pa, 0, &p);
542  geographic_point_init(p.x, p.y, &a);
543 
544  for ( i = 1; i < pa->npoints; i++ )
545  {
546  GEOGRAPHIC_POINT a1, b1;
547  double strip_area = 0.0;
548  double delta_lon = 0.0;
549  LWDEBUGF(4, "edge #%d", i);
550 
551  getPoint2d_p(pa, i, &p);
552  geographic_point_init(p.x, p.y, &b);
553 
554  a1 = a;
555  b1 = b;
556 
557  /* Flip into north if in south */
558  if ( in_south )
559  {
560  a1.lat = -1.0 * a1.lat;
561  b1.lat = -1.0 * b1.lat;
562  }
563 
564  LWDEBUGF(4, "in_south %d", in_south);
565 
566  LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
567 
568  if ( crosses_dateline(&a, &b) )
569  {
570  double shift;
571 
572  if ( a1.lon > 0.0 )
573  shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
574  else
575  shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
576 
577  LWDEBUGF(4, "shift: %.8g", shift);
578  LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
579  point_shift(&a1, shift);
580  point_shift(&b1, shift);
581  LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
582 
583  }
584 
585 
586  delta_lon = fabs(b1.lon - a1.lon);
587 
588  LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
589  LWDEBUGF(4, "delta_lon %.18g", delta_lon);
590  LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
591 
592  if ( delta_lon > 0.0 )
593  {
594  if ( delta_lon < delta_lon_tolerance )
595  {
596  strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
597  LWDEBUGF(4, "strip_area %.12g", strip_area);
598  area += strip_area;
599  }
600  else
601  {
602  GEOGRAPHIC_POINT p, q;
603  double step = floor(delta_lon / delta_lon_tolerance);
604  double distance = spheroid_distance(&a1, &b1, spheroid);
605  double pDistance = 0.0;
606  int j = 0;
607  LWDEBUGF(4, "step %.18g", step);
608  LWDEBUGF(4, "distance %.18g", distance);
609  step = distance / step;
610  LWDEBUGF(4, "step %.18g", step);
611  p = a1;
612  while (pDistance < (distance - step * 1.01))
613  {
614  double azimuth = spheroid_direction(&p, &b1, spheroid);
615  j++;
616  LWDEBUGF(4, " iteration %d", j);
617  LWDEBUGF(4, " azimuth %.12g", azimuth);
618  pDistance = pDistance + step;
619  LWDEBUGF(4, " pDistance %.12g", pDistance);
620  spheroid_project(&p, spheroid, step, azimuth, &q);
621  strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
622  LWDEBUGF(4, " strip_area %.12g", strip_area);
623  area += strip_area;
624  LWDEBUGF(4, " area %.12g", area);
625  p.lat = q.lat;
626  p.lon = q.lon;
627  }
628  strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
629  area += strip_area;
630  }
631  }
632 
633  /* B gets incremented in the next loop, so we save the value here */
634  a = b;
635  }
636  return fabs(area);
637 }
638 #endif /* else ! PROJ_GEODESIC */
639 
647 double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
648 {
649  int type;
650 
651  assert(lwgeom);
652 
653  /* No area in nothing */
654  if ( lwgeom_is_empty(lwgeom) )
655  return 0.0;
656 
657  /* Read the geometry type number */
658  type = lwgeom->type;
659 
660  /* Anything but polygons and collections returns zero */
661  if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
662  return 0.0;
663 
664  /* Actually calculate area */
665  if ( type == POLYGONTYPE )
666  {
667  LWPOLY *poly = (LWPOLY*)lwgeom;
668  uint32_t i;
669  double area = 0.0;
670 
671  /* Just in case there's no rings */
672  if ( poly->nrings < 1 )
673  return 0.0;
674 
675  /* First, the area of the outer ring */
676  area += ptarray_area_spheroid(poly->rings[0], spheroid);
677 
678  /* Subtract areas of inner rings */
679  for ( i = 1; i < poly->nrings; i++ )
680  {
681  area -= ptarray_area_spheroid(poly->rings[i], spheroid);
682  }
683  return area;
684  }
685 
686  /* Recurse into sub-geometries to get area */
687  if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
688  {
689  LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
690  uint32_t i;
691  double area = 0.0;
692 
693  for ( i = 0; i < col->ngeoms; i++ )
694  {
695  area += lwgeom_area_spheroid(col->geoms[i], spheroid);
696  }
697  return area;
698  }
699 
700  /* Shouldn't get here. */
701  return 0.0;
702 }
703 
704 
705 
char * s
Definition: cu_in_wkt.c:23
char * r
Definition: cu_in_wkt.c:24
int ptarray_calculate_gbox_cartesian(const POINTARRAY *pa, GBOX *gbox)
Calculate box (x/y) and add values to gbox.
Definition: g_box.c:542
uint8_t gflags(int hasz, int hasm, int geodetic)
Construct a new flags char.
Definition: g_util.c:145
#define LW_FALSE
Definition: liblwgeom.h:77
#define COLLECTIONTYPE
Definition: liblwgeom.h:91
#define LW_SUCCESS
Definition: liblwgeom.h:80
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:348
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:90
#define POLYGONTYPE
Definition: liblwgeom.h:87
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:1393
#define LW_TRUE
Return types for functions with status returns.
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)
#define FP_MIN(A, B)
void point_shift(GEOGRAPHIC_POINT *p, double shift)
Shift a point around by a number of radians.
Definition: lwgeodetic.c:160
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:180
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:948
int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
Definition: lwgeodetic.c:170
int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
Definition: lwgeodetic.c:666
#define POW2(x)
Definition: lwgeodetic.h:47
#define deg2rad(d)
Conversion functions.
Definition: lwgeodetic.h:79
Datum distance(PG_FUNCTION_ARGS)
Datum area(PG_FUNCTION_ARGS)
#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:190
double deltaLongitude(double azimuth, double sigma, double tsm, SPHEROID *sphere)
static double spheroid_big_b(double u2)
Definition: lwspheroid.c:60
static double spheroid_big_a(double u2)
Definition: lwspheroid.c:55
static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
Definition: lwspheroid.c:499
static double spheroid_mu2(double alpha, const SPHEROID *s)
Definition: lwspheroid.c:49
double lwgeom_area_spheroid(const LWGEOM *lwgeom, const SPHEROID *spheroid)
Calculate the area of an LWGEOM.
Definition: lwspheroid.c:647
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:191
static double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
Definition: lwspheroid.c:425
void spheroid_init(SPHEROID *s, double a, double b)
Initialize spheroid object based on major and minor axis.
Definition: lwspheroid.c:39
static double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
Definition: lwspheroid.c:420
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:460
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:442
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:364
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:287
type
Definition: ovdump.py:41
double ymax
Definition: liblwgeom.h:298
double ymin
Definition: liblwgeom.h:297
uint8_t flags
Definition: liblwgeom.h:294
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:53
uint32_t ngeoms
Definition: liblwgeom.h:510
LWGEOM ** geoms
Definition: liblwgeom.h:512
uint8_t type
Definition: liblwgeom.h:399
POINTARRAY ** rings
Definition: liblwgeom.h:460
uint32_t nrings
Definition: liblwgeom.h:458
double y
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:331
uint32_t npoints
Definition: liblwgeom.h:374
double e_sq
Definition: liblwgeom.h:320
double radius
Definition: liblwgeom.h:321
double a
Definition: liblwgeom.h:316
double b
Definition: liblwgeom.h:317
double f
Definition: liblwgeom.h:318
unsigned int uint32_t
Definition: uthash.h:78