00001 #include <grass/gis.h>
00002
00003 static struct Cell_head window;
00004 static double square_meters;
00005 static int projection;
00006
00007 static double units_to_meters_squared = 0.0;
00008
00009
00010 static int next_row;
00011 static double north_value;
00012 static double north;
00013 static double (*darea0)();
00014
00015
00033 int G_begin_cell_area_calculations()
00034 {
00035 double a, e2;
00036 double factor;
00037
00038 G_get_set_window(&window);
00039 switch(projection = window.proj)
00040 {
00041 case PROJECTION_LL:
00042 G_get_ellipsoid_parameters (&a, &e2);
00043 if (e2)
00044 {
00045 G_begin_zone_area_on_ellipsoid (a, e2, window.ew_res/360.0);
00046 darea0 = G_darea0_on_ellipsoid;
00047 }
00048 else
00049 {
00050 G_begin_zone_area_on_sphere (a, window.ew_res/360.0);
00051 darea0 = G_darea0_on_sphere;
00052 }
00053 next_row = 0;
00054 north_value = darea0 (north = window.north);
00055 return 2;
00056 default:
00057 square_meters = window.ns_res * window.ew_res;
00058 factor = G_database_units_to_meters_factor();
00059 if (factor > 0.0)
00060 square_meters *= (factor * factor);
00061 return (factor > 0.0);
00062 }
00063 }
00064
00065
00077 double
00078 G_area_of_cell_at_row ( register int row)
00079 {
00080 register double south_value;
00081 register double cell_area;
00082
00083 if (projection != PROJECTION_LL)
00084 return square_meters;
00085
00086 if (row != next_row)
00087 north_value = darea0 (north = window.north - row * window.ns_res);
00088
00089 south_value = darea0 (north -= window.ns_res);
00090 cell_area = north_value - south_value;
00091
00092 next_row = row+1;
00093 north_value = south_value;
00094
00095 return cell_area;
00096 }
00097
00098
00113 int G_begin_polygon_area_calculations()
00114 {
00115 double a, e2;
00116 double factor;
00117
00118 if ((projection = G_projection()) == PROJECTION_LL)
00119 {
00120 G_get_ellipsoid_parameters (&a, &e2);
00121 G_begin_ellipsoid_polygon_area (a, e2);
00122 return 2;
00123 }
00124 factor = G_database_units_to_meters_factor();
00125 if (factor > 0.0)
00126 {
00127 units_to_meters_squared = factor *factor;
00128 return 1;
00129 }
00130 units_to_meters_squared = 1.0;
00131 return 0;
00132 }
00133
00134
00156 double G_area_of_polygon(double *x,double *y,int n)
00157 {
00158 double area;
00159
00160 if (projection == PROJECTION_LL)
00161 area = G_ellipsoid_polygon_area(x,y,n);
00162 else
00163 area = G_planimetric_polygon_area(x,y,n) * units_to_meters_squared;
00164
00165 return area;
00166 }