source: trunk/liblwgeom/lwspheroid.c @ 5274

Last change on this file since 5274 was 5274, checked in by colivier, 6 years ago

make astyle session

  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1/**********************************************************************
2 * $Id: lwspheroid.c 5274 2010-02-21 12:36:27Z colivier $
3 *
4 * PostGIS - Spatial Types for PostgreSQL
5 * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
6 * Copyright 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 "lwgeodetic.h"
14
15/**
16* Initialize spheroid object based on major and minor axis
17*/
18void spheroid_init(SPHEROID *s, double a, double b)
19{
20        s->a = a;
21        s->b = b;
22        s->f = (a - b) / a;
23        s->e_sq = (a*a - b*b)/(a*a);
24        s->radius = (2.0 * a + b ) / 3.0;
25}
26
27static double spheroid_mu2(double alpha, const SPHEROID *s)
28{
29        double b2 = POW2(s->b);
30        return POW2(cos(alpha)) * (POW2(s->a) - b2) / b2;
31}
32
33static double spheroid_big_a(double u2)
34{
35        return 1.0 + (u2 / 16384.0) * (4096.0 + u2 * (-768.0 + u2 * (320.0 - 175.0 * u2)));
36}
37
38static double spheroid_big_b(double u2)
39{
40        return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
41}
42
43/**
44* Computes the shortest distance along the surface of the spheroid
45* between two points. Based on Vincenty's formula for the geodetic
46* inverse problem as described in "Geocentric Datum of Australia
47* Technical Manual", Chapter 4. Tested against:
48* http://mascot.gdbc.gov.bc.ca/mascot/util1a.html
49* and
50* http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp
51*
52* @param a - location of first point.
53* @param b - location of second point.
54* @param s - spheroid to calculate on
55* @return spheroidal distance between a and b in spheroid units.
56*/
57double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
58{
59        double lambda = (b->lon - a->lon);
60        double f = spheroid->f;
61        double omf = 1 - spheroid->f;
62        double u1, u2;
63        double cos_u1, cos_u2;
64        double sin_u1, sin_u2;
65        double big_a, big_b, delta_sigma;
66        double alpha, sin_alpha, cos_alphasq, c;
67        double sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqrsin_sigma, last_lambda, omega;
68        double cos_lambda, sin_lambda;
69        double distance;
70        int i = 0;
71
72        /* Same point => zero distance */
73        if ( geographic_point_equals(a, b) )
74        {
75                return 0.0;
76        }
77
78        u1 = atan(omf * tan(a->lat));
79        cos_u1 = cos(u1);
80        sin_u1 = sin(u1);
81        u2 = atan(omf * tan(b->lat));
82        cos_u2 = cos(u2);
83        sin_u2 = sin(u2);
84
85        omega = lambda;
86        do
87        {
88                cos_lambda = cos(lambda);
89                sin_lambda = sin(lambda);
90                sqrsin_sigma = POW2(cos_u2 * sin_lambda) +
91                               POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda));
92                sin_sigma = sqrt(sqrsin_sigma);
93                cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
94                sigma = atan2(sin_sigma, cos_sigma);
95                sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin(sigma);
96
97                /* Numerical stability issue, ensure asin is not NaN */
98                if ( sin_alpha > 1.0 )
99                        alpha = M_PI_2;
100                else if ( sin_alpha < -1.0 )
101                        alpha = -1.0 * M_PI_2;
102                else
103                        alpha = asin(sin_alpha);
104
105                cos_alphasq = POW2(cos(alpha));
106                cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
107
108                /* Numerical stability issue, cos2 is in range */
109                if ( cos2_sigma_m > 1.0 )
110                        cos2_sigma_m = 1.0;
111                if ( cos2_sigma_m < -1.0 )
112                        cos2_sigma_m = -1.0;
113
114                c = (f / 16.0) * cos_alphasq * (4.0 + f * (4.0 - 3.0 * cos_alphasq));
115                last_lambda = lambda;
116                lambda = omega + (1.0 - c) * f * sin(alpha) * (sigma + c * sin(sigma) *
117                         (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
118                i++;
119        }
120        while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
121
122        u2 = spheroid_mu2(alpha, spheroid);
123        big_a = spheroid_big_a(u2);
124        big_b = spheroid_big_b(u2);
125        delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
126                                           (big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
127
128        distance = spheroid->b * big_a * (sigma - delta_sigma);
129
130        /* Algorithm failure, distance == NaN, fallback to sphere */
131        if ( distance != distance )
132        {
133                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);
134                return spheroid->radius * sphere_distance(a, b);
135        }
136
137        return distance;
138}
139
140/**
141* Computes the direction of the geodesic joining two points on
142* the spheroid. Based on Vincenty's formula for the geodetic
143* inverse problem as described in "Geocentric Datum of Australia
144* Technical Manual", Chapter 4. Tested against:
145* http://mascot.gdbc.gov.bc.ca/mascot/util1a.html
146* and
147* http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp
148*
149* @param r - location of first point
150* @param s - location of second point
151* @return azimuth of line joining r and s
152*/
153double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
154{
155        int i = 0;
156        double lambda = s->lon - r->lon;
157        double omf = 1 - spheroid->f;
158        double u1 = atan(omf * tan(r->lat));
159        double cos_u1 = cos(u1);
160        double sin_u1 = sin(u1);
161        double u2 = atan(omf * tan(s->lat));
162        double cos_u2 = cos(u2);
163        double sin_u2 = sin(u2);
164
165        double omega = lambda;
166        double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
167        double sin_alpha, cos_alphasq, C, alphaFD;
168        do
169        {
170                sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
171                                POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
172                sin_sigma = sqrt(sqr_sin_sigma);
173                cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
174                sigma = atan2(sin_sigma, cos_sigma);
175                sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
176
177                /* Numerical stability issue, ensure asin is not NaN */
178                if ( sin_alpha > 1.0 )
179                        alpha = M_PI_2;
180                else if ( sin_alpha < -1.0 )
181                        alpha = -1.0 * M_PI_2;
182                else
183                        alpha = asin(sin_alpha);
184
185                alpha = asin(sin_alpha);
186                cos_alphasq = POW2(cos(alpha));
187                cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
188
189                /* Numerical stability issue, cos2 is in range */
190                if ( cos2_sigma_m > 1.0 )
191                        cos2_sigma_m = 1.0;
192                if ( cos2_sigma_m < -1.0 )
193                        cos2_sigma_m = -1.0;
194
195                C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
196                last_lambda = lambda;
197                lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
198                         (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
199                i++;
200        }
201        while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
202
203        alphaFD = atan2((cos_u2 * sin(lambda)),
204                        (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
205        if (alphaFD < 0.0)
206        {
207                alphaFD = alphaFD + 2.0 * M_PI;
208        }
209        if (alphaFD > 2.0 * M_PI)
210        {
211                alphaFD = alphaFD - 2.0 * M_PI;
212        }
213        return alphaFD;
214}
215
216
217/**
218* Given a location, an azimuth and a distance, computes the
219* location of the projected point. Based on Vincenty's formula
220* for the geodetic direct problem as described in "Geocentric
221* Datum of Australia Technical Manual", Chapter 4. Tested against:
222* http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
223* and
224* http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
225*
226* @param r - location of first point.
227* @param distance - distance in meters.
228* @param azimuth - azimuth in radians.
229* @return s - location of projected point.
230*/
231int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
232{
233        double omf = 1 - spheroid->f;
234        double tan_u1 = omf * tan(r->lat);
235        double u1 = atan(tan_u1);
236        double sigma, last_sigma, delta_sigma, two_sigma_m;
237        double sigma1, sin_alpha, alpha, cos_alphasq;
238        double u2, A, B;
239        double lat2, lambda, lambda2, C, omega;
240        int i = 0;
241
242        if (azimuth < 0.0)
243        {
244                azimuth = azimuth + M_PI * 2.0;
245        }
246        if (azimuth > (PI * 2.0))
247        {
248                azimuth = azimuth - M_PI * 2.0;
249        }
250
251        sigma1 = atan2(tan_u1, cos(azimuth));
252        sin_alpha = cos(u1) * sin(azimuth);
253        alpha = asin(sin_alpha);
254        cos_alphasq = 1.0 - POW2(sin_alpha);
255
256        u2 = spheroid_mu2(alpha, spheroid);
257        A = spheroid_big_a(u2);
258        B = spheroid_big_b(u2);
259
260        sigma = (distance / (spheroid->b * A));
261        do
262        {
263                two_sigma_m = 2.0 * sigma1 + sigma;
264                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))))));
265                last_sigma = sigma;
266                sigma = (distance / (spheroid->b * A)) + delta_sigma;
267                i++;
268        }
269        while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
270
271        lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
272                      cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
273                                      POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
274                                           cos(azimuth)))));
275        lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
276                       sin(u1) * sin(sigma) * cos(azimuth)));
277        C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
278        omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
279                (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
280        lambda2 = r->lon + omega;
281        g->lat = lat2;
282        g->lon = lambda2;
283        return G_SUCCESS;
284}
285
286
287static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
288{
289        return spheroid->a / (sqrt(1.0 - spheroid->e_sq * POW2(sin(latitude))));
290}
291
292static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
293{
294        return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid)
295               * cos(latitude)
296               * deltaLongitude;
297}
298
299
300/**
301* Computes the area on the spheroid of a box bounded by meridians and
302* parallels. The box is defined by two points, the South West corner
303* and the North East corner. Formula based on Bagratuni 1967.
304*
305* @param southWestCorner - lower left corner of bounding box.
306* @param northEastCorner - upper right corner of bounding box.
307* @return area in square meters.
308*/
309static double spheroid_boundingbox_area(const GEOGRAPHIC_POINT *southWestCorner, const GEOGRAPHIC_POINT *northEastCorner, const SPHEROID *spheroid)
310{
311        double z0 = (northEastCorner->lon - southWestCorner->lon) * POW2(spheroid->b) / 2.0;
312        double e = sqrt(spheroid->e_sq);
313        double sinPhi1 = sin(southWestCorner->lat);
314        double sinPhi2 = sin(northEastCorner->lat);
315        double t1p1 = sinPhi1 / (1.0 - spheroid->e_sq * sinPhi1 * sinPhi1);
316        double t1p2 = sinPhi2 / (1.0 - spheroid->e_sq * sinPhi2 * sinPhi2);
317        double oneOver2e = 1.0 / (2.0 * e);
318        double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1));
319        double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2));
320        return z0 * (t1p2 + t2p2) - z0 * (t1p1 + t2p1);
321}
322
323/**
324* This function doesn't work for edges crossing the dateline or in the southern
325* hemisphere. Points are pre-conditioned in ptarray_area_spheroid.
326*/
327static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
328{
329        GEOGRAPHIC_POINT A, B, mL, nR;
330        double deltaLng, baseArea, topArea;
331        double bE, tE, ratio, sign;
332
333        A = *a;
334        B = *b;
335
336        mL.lat = latitude_min;
337        mL.lon = FP_MIN(A.lon, B.lon);
338        nR.lat = FP_MIN(A.lat, B.lat);
339        nR.lon = FP_MAX(A.lon, B.lon);
340        LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
341        LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
342        baseArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
343        LWDEBUGF(4, "baseArea %.12g", baseArea);
344
345        mL.lat = FP_MIN(A.lat, B.lat);
346        mL.lon = FP_MIN(A.lon, B.lon);
347        nR.lat = FP_MAX(A.lat, B.lat);
348        nR.lon = FP_MAX(A.lon, B.lon);
349        LWDEBUGF(4, "mL (%.12g %.12g)", mL.lat, mL.lon);
350        LWDEBUGF(4, "nR (%.12g %.12g)", nR.lat, nR.lon);
351        topArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
352        LWDEBUGF(4, "topArea %.12g", topArea);
353
354        deltaLng = B.lon - A.lon;
355        LWDEBUGF(4, "deltaLng %.12g", deltaLng);
356        bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid);
357        tE = spheroid_parallel_arc_length(B.lat, deltaLng, spheroid);
358        LWDEBUGF(4, "bE %.12g", bE);
359        LWDEBUGF(4, "tE %.12g", tE);
360
361        ratio = (bE + tE)/tE;
362        sign = signum(B.lon - A.lon);
363        return (baseArea + topArea / ratio) * sign;
364}
365
366static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
367{
368        GEOGRAPHIC_POINT a, b;
369        POINT2D p;
370        int i;
371        double area = 0.0;
372        GBOX gbox2d;
373        int in_south = LW_FALSE;
374        double delta_lon_tolerance;
375        double latitude_min;
376
377        gbox2d.flags = gflags(0, 0, 0);
378
379        /* Return zero on non-sensical inputs */
380        if ( ! pa || pa->npoints < 4 )
381                return 0.0;
382
383        /* Get the raw min/max values for the latitudes */
384        ptarray_calculate_gbox(pa, &gbox2d);
385
386        if ( signum(gbox2d.ymin) != signum(gbox2d.ymax) )
387                lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
388
389        /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
390        if ( gbox2d.ymax < 0.0 )
391                in_south = LW_TRUE;
392
393        LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
394
395        /* Tolerance for strip area calculation */
396        if ( in_south )
397        {
398                delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
399                latitude_min = deg2rad(fabs(gbox2d.ymax));
400        }
401        else
402        {
403                delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
404                latitude_min = deg2rad(gbox2d.ymin);
405        }
406
407        /* Initialize first point */
408        getPoint2d_p(pa, 0, &p);
409        geographic_point_init(p.x, p.y, &a);
410
411        for ( i = 1; i < pa->npoints; i++ )
412        {
413                GEOGRAPHIC_POINT a1, b1;
414                double strip_area = 0.0;
415                double delta_lon = 0.0;
416                LWDEBUGF(4, "edge #%d", i);
417
418                getPoint2d_p(pa, i, &p);
419                geographic_point_init(p.x, p.y, &b);
420
421                a1 = a;
422                b1 = b;
423
424                /* Flip into north if in south */
425                if ( in_south )
426                {
427                        a1.lat = -1.0 * a1.lat;
428                        b1.lat = -1.0 * b1.lat;
429                }
430
431                LWDEBUGF(4, "in_south %d", in_south);
432
433                if ( crosses_dateline(&a, &b) )
434                {
435                        double shift;
436
437                        if ( a1.lon > 0.0 )
438                                shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
439                        else
440                                shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
441
442                        point_shift(&a1, shift);
443                        point_shift(&b1, shift);
444                }
445
446                LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
447
448                delta_lon = fabs(b1.lon - a1.lon);
449
450                LWDEBUGF(4, "(%.18g %.18g) (%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
451                LWDEBUGF(4, "delta_lon %.18g", delta_lon);
452                LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
453
454                if ( delta_lon > 0.0 )
455                {
456                        if ( delta_lon < delta_lon_tolerance )
457                        {
458                                strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
459                                LWDEBUGF(4, "strip_area %.12g", strip_area);
460                                area += strip_area;
461                        }
462                        else
463                        {
464                                GEOGRAPHIC_POINT p, q;
465                                double step = floor(delta_lon / delta_lon_tolerance);
466                                double distance = spheroid_distance(&a1, &b1, spheroid);
467                                double pDistance = 0.0;
468                                int j = 0;
469                                LWDEBUGF(4, "step %.18g", step);
470                                LWDEBUGF(4, "distance %.18g", distance);
471                                step = distance / step;
472                                LWDEBUGF(4, "step %.18g", step);
473                                p = a1;
474                                while (pDistance < (distance - step * 1.01))
475                                {
476                                        double azimuth = spheroid_direction(&p, &b1, spheroid);
477                                        j++;
478                                        LWDEBUGF(4, "  iteration %d", j);
479                                        LWDEBUGF(4, "  azimuth %.12g", azimuth);
480                                        pDistance = pDistance + step;
481                                        LWDEBUGF(4, "  pDistance %.12g", pDistance);
482                                        spheroid_project(&p, spheroid, step, azimuth, &q);
483                                        strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
484                                        LWDEBUGF(4, "  strip_area %.12g", strip_area);
485                                        area += strip_area;
486                                        LWDEBUGF(4, "  area %.12g", area);
487                                        p.lat = q.lat;
488                                        p.lon = q.lon;
489                                }
490                                strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
491                                area += strip_area;
492                        }
493                }
494
495                /* B gets incremented in the next loop, so we save the value here */
496                a = b;
497        }
498        return fabs(area);
499}
500/**
501* Calculate the area of an LWGEOM. Anything except POLYGON, MULTIPOLYGON
502* and GEOMETRYCOLLECTION return zero immediately. Multi's recurse, polygons
503* calculate external ring area and subtract internal ring area. A GBOX is
504* required to check relationship to equator an outside point.
505* WARNING: Does NOT WORK for polygons over equator or pole.
506*/
507double lwgeom_area_spheroid(const LWGEOM *lwgeom, const GBOX *gbox, const SPHEROID *spheroid)
508{
509        int type;
510
511        assert(lwgeom);
512
513        /* No area in nothing */
514        if ( lwgeom_is_empty(lwgeom) )
515                return 0.0;
516
517        /* Read the geometry type number */
518        type = TYPE_GETTYPE(lwgeom->type);
519
520        /* Anything but polygons and collections returns zero */
521        if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
522                return 0.0;
523
524        /* Actually calculate area */
525        if ( type == POLYGONTYPE )
526        {
527                LWPOLY *poly = (LWPOLY*)lwgeom;
528                int i;
529                double area = 0.0;
530
531                /* Just in case there's no rings */
532                if ( poly->nrings < 1 )
533                        return 0.0;
534
535                /* First, the area of the outer ring */
536                area += ptarray_area_spheroid(poly->rings[0], spheroid);
537
538                /* Subtract areas of inner rings */
539                for ( i = 1; i < poly->nrings; i++ )
540                {
541                        area -=  ptarray_area_spheroid(poly->rings[i], spheroid);
542                }
543                return area;
544        }
545
546        /* Recurse into sub-geometries to get area */
547        if ( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
548        {
549                LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
550                int i;
551                double area = 0.0;
552
553                for ( i = 0; i < col->ngeoms; i++ )
554                {
555                        area += lwgeom_area_spheroid(col->geoms[i], gbox, spheroid);
556                }
557                return area;
558        }
559
560        /* Shouldn't get here. */
561        return 0.0;
562}
563
564
565
Note: See TracBrowser for help on using the repository browser.