ptarray_area_spheroid()

static double ptarray_area_spheroid ( const POINTARRAY *  pa, const SPHEROID *  spheroid  )

static

Definition at line 494 of file lwspheroid.c.

References area(), crosses_dateline(), deg2rad, distance(), GBOX::flags, geographic_point_init(), getPoint2d_p(), gflags(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, LW_FALSE, LW_TRUE, LWDEBUGF, lwerror(), POINTARRAY::npoints, point_shift(), ptarray_calculate_gbox_cartesian(), SIGNUM, spheroid_direction(), spheroid_distance(), spheroid_project(), spheroid_striparea(), POINT2D::x, POINT2D::y, GBOX::ymax, and GBOX::ymin.

Referenced by lwgeom_area_spheroid(), and spheroid_big_b().

{
        GEOGRAPHIC_POINT a, b;
        POINT2D p;
        int i;
        double area = 0.0;
        GBOX gbox2d;
        int in_south = LW_FALSE;
        double delta_lon_tolerance;
        double latitude_min;

        gbox2d.flags = gflags(0, 0, 0);

        /* Return zero on non-sensical inputs */
        if ( ! pa || pa->npoints < 4 )
                return 0.0;

        /* Get the raw min/max values for the latitudes */
        ptarray_calculate_gbox_cartesian(pa, &gbox2d);

        if ( SIGNUM(gbox2d.ymin) != SIGNUM(gbox2d.ymax) )
                lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");

        /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
        if ( gbox2d.ymax < 0.0 )
                in_south = LW_TRUE;

        LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);

        /* Tolerance for strip area calculation */
        if ( in_south )
        {
                delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
                latitude_min = deg2rad(fabs(gbox2d.ymax));
        }
        else
        {
                delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
                latitude_min = deg2rad(gbox2d.ymin);
        }

        /* Initialize first point */
        getPoint2d_p(pa, 0, &p);
        geographic_point_init(p.x, p.y, &a);

        for ( i = 1; i < pa->npoints; i++ )
        {
                GEOGRAPHIC_POINT a1, b1;
                double strip_area = 0.0;
                double delta_lon = 0.0;
                LWDEBUGF(4, "edge #%d", i);

                getPoint2d_p(pa, i, &p);
                geographic_point_init(p.x, p.y, &b);

                a1 = a;
                b1 = b;

                /* Flip into north if in south */
                if ( in_south )
                {
                        a1.lat = -1.0 * a1.lat;
                        b1.lat = -1.0 * b1.lat;
                }

                LWDEBUGF(4, "in_south %d", in_south);

                LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );

                if ( crosses_dateline(&a, &b) )
                {
                        double shift;

                        if ( a1.lon > 0.0 )
                                shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
                        else
                                shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */

                        LWDEBUGF(4, "shift: %.8g", shift);
                        LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
                        point_shift(&a1, shift);
                        point_shift(&b1, shift);
                        LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);

                }


                delta_lon = fabs(b1.lon - a1.lon);

                LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
                LWDEBUGF(4, "delta_lon %.18g", delta_lon);
                LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);

                if ( delta_lon > 0.0 )
                {
                        if ( delta_lon < delta_lon_tolerance )
                        {
                                strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
                                LWDEBUGF(4, "strip_area %.12g", strip_area);
                                area += strip_area;
                        }
                        else
                        {
                                GEOGRAPHIC_POINT p, q;
                                double step = floor(delta_lon / delta_lon_tolerance);
                                double distance = spheroid_distance(&a1, &b1, spheroid);
                                double pDistance = 0.0;
                                int j = 0;
                                LWDEBUGF(4, "step %.18g", step);
                                LWDEBUGF(4, "distance %.18g", distance);
                                step = distance / step;
                                LWDEBUGF(4, "step %.18g", step);
                                p = a1;
                                while (pDistance < (distance - step * 1.01))
                                {
                                        double azimuth = spheroid_direction(&p, &b1, spheroid);
                                        j++;
                                        LWDEBUGF(4, "  iteration %d", j);
                                        LWDEBUGF(4, "  azimuth %.12g", azimuth);
                                        pDistance = pDistance + step;
                                        LWDEBUGF(4, "  pDistance %.12g", pDistance);
                                        spheroid_project(&p, spheroid, step, azimuth, &q);
                                        strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
                                        LWDEBUGF(4, "  strip_area %.12g", strip_area);
                                        area += strip_area;
                                        LWDEBUGF(4, "  area %.12g", area);
                                        p.lat = q.lat;
                                        p.lon = q.lon;
                                }
                                strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
                                area += strip_area;
                        }
                }

                /* B gets incremented in the next loop, so we save the value here */
                a = b;
        }
        return fabs(area);
}

Here is the call graph for this function:

Here is the caller graph for this function: