00001
00017 #include <grass/gis.h>
00018
00019
00020 static struct Cell_head window;
00021 static double square_meters;
00022 static int projection;
00023
00024 static double units_to_meters_squared = 0.0;
00025
00026
00027 static int next_row;
00028 static double north_value;
00029 static double north;
00030 static double (*darea0) (double);
00031
00032
00050 int G_begin_cell_area_calculations(void)
00051 {
00052 double a, e2;
00053 double factor;
00054
00055 G_get_set_window(&window);
00056 switch (projection = window.proj) {
00057 case PROJECTION_LL:
00058 G_get_ellipsoid_parameters(&a, &e2);
00059 if (e2) {
00060 G_begin_zone_area_on_ellipsoid(a, e2, window.ew_res / 360.0);
00061 darea0 = G_darea0_on_ellipsoid;
00062 }
00063 else {
00064 G_begin_zone_area_on_sphere(a, window.ew_res / 360.0);
00065 darea0 = G_darea0_on_sphere;
00066 }
00067 next_row = 0;
00068 north_value = darea0(north = window.north);
00069 return 2;
00070 default:
00071 square_meters = window.ns_res * window.ew_res;
00072 factor = G_database_units_to_meters_factor();
00073 if (factor > 0.0)
00074 square_meters *= (factor * factor);
00075 return (factor > 0.0);
00076 }
00077 }
00078
00079
00091 double G_area_of_cell_at_row(int row)
00092 {
00093 register double south_value;
00094 register double cell_area;
00095
00096 if (projection != PROJECTION_LL)
00097 return square_meters;
00098
00099 if (row != next_row)
00100 north_value = darea0(north = window.north - row * window.ns_res);
00101
00102 south_value = darea0(north -= window.ns_res);
00103 cell_area = north_value - south_value;
00104
00105 next_row = row + 1;
00106 north_value = south_value;
00107
00108 return cell_area;
00109 }
00110
00111
00123 int G_begin_polygon_area_calculations(void)
00124 {
00125 double a, e2;
00126 double factor;
00127
00128 if ((projection = G_projection()) == PROJECTION_LL) {
00129 G_get_ellipsoid_parameters(&a, &e2);
00130 G_begin_ellipsoid_polygon_area(a, e2);
00131 return 2;
00132 }
00133 factor = G_database_units_to_meters_factor();
00134 if (factor > 0.0) {
00135 units_to_meters_squared = factor * factor;
00136 return 1;
00137 }
00138 units_to_meters_squared = 1.0;
00139 return 0;
00140 }
00141
00142
00165 double G_area_of_polygon(const double *x, const double *y, int n)
00166 {
00167 double area;
00168
00169 if (projection == PROJECTION_LL)
00170 area = G_ellipsoid_polygon_area(x, y, n);
00171 else
00172 area = G_planimetric_polygon_area(x, y, n) * units_to_meters_squared;
00173
00174 return area;
00175 }