static double ptarray_area_spheroid ( const POINTARRAY pa,
const SPHEROID spheroid 
) [static]

Definition at line 367 of file lwspheroid.c.

References a, crosses_dateline(), deg2rad, 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().

00368 {
00369         GEOGRAPHIC_POINT a, b;
00370         POINT2D p;
00371         int i;
00372         double area = 0.0;
00373         GBOX gbox2d;
00374         int in_south = LW_FALSE;
00375         double delta_lon_tolerance;
00376         double latitude_min;
00377 
00378         gbox2d.flags = gflags(0, 0, 0);
00379 
00380         /* Return zero on non-sensical inputs */
00381         if ( ! pa || pa->npoints < 4 )
00382                 return 0.0;
00383 
00384         /* Get the raw min/max values for the latitudes */
00385         ptarray_calculate_gbox_cartesian(pa, &gbox2d);
00386 
00387         if ( signum(gbox2d.ymin) != signum(gbox2d.ymax) )
00388                 lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
00389 
00390         /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
00391         if ( gbox2d.ymax < 0.0 )
00392                 in_south = LW_TRUE;
00393 
00394         LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
00395 
00396         /* Tolerance for strip area calculation */
00397         if ( in_south )
00398         {
00399                 delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
00400                 latitude_min = deg2rad(fabs(gbox2d.ymax));
00401         }
00402         else
00403         {
00404                 delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
00405                 latitude_min = deg2rad(gbox2d.ymin);
00406         }
00407 
00408         /* Initialize first point */
00409         getPoint2d_p(pa, 0, &p);
00410         geographic_point_init(p.x, p.y, &a);
00411 
00412         for ( i = 1; i < pa->npoints; i++ )
00413         {
00414                 GEOGRAPHIC_POINT a1, b1;
00415                 double strip_area = 0.0;
00416                 double delta_lon = 0.0;
00417                 LWDEBUGF(4, "edge #%d", i);
00418 
00419                 getPoint2d_p(pa, i, &p);
00420                 geographic_point_init(p.x, p.y, &b);
00421 
00422                 a1 = a;
00423                 b1 = b;
00424 
00425                 /* Flip into north if in south */
00426                 if ( in_south )
00427                 {
00428                         a1.lat = -1.0 * a1.lat;
00429                         b1.lat = -1.0 * b1.lat;
00430                 }
00431 
00432                 LWDEBUGF(4, "in_south %d", in_south);
00433 
00434                 LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
00435 
00436                 if ( crosses_dateline(&a, &b) )
00437                 {
00438                         double shift;
00439 
00440                         if ( a1.lon > 0.0 )
00441                                 shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
00442                         else
00443                                 shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
00444 
00445                         LWDEBUGF(4, "shift: %.8g", shift);
00446                         LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
00447                         point_shift(&a1, shift);
00448                         point_shift(&b1, shift);
00449                         LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon);
00450 
00451                 }
00452 
00453 
00454                 delta_lon = fabs(b1.lon - a1.lon);
00455 
00456                 LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
00457                 LWDEBUGF(4, "delta_lon %.18g", delta_lon);
00458                 LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
00459 
00460                 if ( delta_lon > 0.0 )
00461                 {
00462                         if ( delta_lon < delta_lon_tolerance )
00463                         {
00464                                 strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
00465                                 LWDEBUGF(4, "strip_area %.12g", strip_area);
00466                                 area += strip_area;
00467                         }
00468                         else
00469                         {
00470                                 GEOGRAPHIC_POINT p, q;
00471                                 double step = floor(delta_lon / delta_lon_tolerance);
00472                                 double distance = spheroid_distance(&a1, &b1, spheroid);
00473                                 double pDistance = 0.0;
00474                                 int j = 0;
00475                                 LWDEBUGF(4, "step %.18g", step);
00476                                 LWDEBUGF(4, "distance %.18g", distance);
00477                                 step = distance / step;
00478                                 LWDEBUGF(4, "step %.18g", step);
00479                                 p = a1;
00480                                 while (pDistance < (distance - step * 1.01))
00481                                 {
00482                                         double azimuth = spheroid_direction(&p, &b1, spheroid);
00483                                         j++;
00484                                         LWDEBUGF(4, "  iteration %d", j);
00485                                         LWDEBUGF(4, "  azimuth %.12g", azimuth);
00486                                         pDistance = pDistance + step;
00487                                         LWDEBUGF(4, "  pDistance %.12g", pDistance);
00488                                         spheroid_project(&p, spheroid, step, azimuth, &q);
00489                                         strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
00490                                         LWDEBUGF(4, "  strip_area %.12g", strip_area);
00491                                         area += strip_area;
00492                                         LWDEBUGF(4, "  area %.12g", area);
00493                                         p.lat = q.lat;
00494                                         p.lon = q.lon;
00495                                 }
00496                                 strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
00497                                 area += strip_area;
00498                         }
00499                 }
00500 
00501                 /* B gets incremented in the next loop, so we save the value here */
00502                 a = b;
00503         }
00504         return fabs(area);
00505 }

Here is the call graph for this function:

Here is the caller graph for this function:


Generated on Wed Aug 22 01:36:30 2012 for PostGIS Trunk Doxygen by  doxygen 1.4.7