GRASS Programmer's Manual
6.4.2(2012)
|
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 /* these next are for lat-long only */ 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 }