GRASS Programmer's Manual
6.4.1(2011)
|
00001 00017 #include <grass/gis.h> 00018 #include <grass/glocale.h> 00019 00020 00043 char *G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag) 00044 { 00045 if (!row_flag) { 00046 if (cellhd->ns_res <= 0) 00047 return (_("Illegal n-s resolution value")); 00048 } 00049 else { 00050 if (cellhd->rows <= 0) 00051 return (_("Illegal row value")); 00052 } 00053 if (!col_flag) { 00054 if (cellhd->ew_res <= 0) 00055 return (_("Illegal e-w resolution value")); 00056 } 00057 else { 00058 if (cellhd->cols <= 0) 00059 return (_("Illegal col value")); 00060 } 00061 00062 /* for lat/lon, check north,south. force east larger than west */ 00063 if (cellhd->proj == PROJECTION_LL) { 00064 double epsilon_ns, epsilon_ew; 00065 00066 /* TODO: find good thresholds */ 00067 epsilon_ns = 1. / cellhd->rows * 0.001; 00068 epsilon_ew = .000001; /* epsilon_ew calculation doesn't work due to cellhd->cols update/global wraparound below */ 00069 00070 G_debug(3, "G_adjust_Cell_head: epsilon_ns: %g, epsilon_ew: %g", 00071 epsilon_ns, epsilon_ew); 00072 00073 /* TODO: once working, change below G_warning to G_debug */ 00074 00075 /* fix rounding problems if input map slightly exceeds the world definition -180 90 180 -90 */ 00076 if (cellhd->north > 90.0) { 00077 if (((cellhd->north - 90.0) < epsilon_ns) && 00078 ((cellhd->north - 90.0) > GRASS_EPSILON)) { 00079 G_warning(_("Fixing subtle input data rounding error of north boundary (%g>%g)"), 00080 cellhd->north - 90.0, epsilon_ns); 00081 cellhd->north = 90.0; 00082 } 00083 else 00084 return (_("Illegal latitude for North")); 00085 } 00086 00087 if (cellhd->south < -90.0) { 00088 if (((cellhd->south + 90.0) < epsilon_ns) && 00089 ((cellhd->south + 90.0) < GRASS_EPSILON)) { 00090 G_warning(_("Fixing subtle input data rounding error of south boundary (%g>%g)"), 00091 cellhd->south + 90.0, epsilon_ns); 00092 cellhd->south = -90.0; 00093 } 00094 else 00095 return (_("Illegal latitude for South")); 00096 } 00097 00098 #if 0 00099 /* DISABLED: this breaks global wrap-around */ 00100 00101 G_debug(3, 00102 "G_adjust_Cell_head() cellhd->west: %f, devi: %g, eps: %g", 00103 cellhd->west, cellhd->west + 180.0, epsilon_ew); 00104 00105 if ((cellhd->west < -180.0) && ((cellhd->west + 180.0) < epsilon_ew) 00106 && ((cellhd->west + 180.0) < GRASS_EPSILON)) { 00107 G_warning(_("Fixing subtle input data rounding error of west boundary (%g>%g)"), 00108 cellhd->west + 180.0, epsilon_ew); 00109 cellhd->west = -180.0; 00110 } 00111 00112 G_debug(3, 00113 "G_adjust_Cell_head() cellhd->east: %f, devi: %g, eps: %g", 00114 cellhd->east, cellhd->east - 180.0, epsilon_ew); 00115 00116 if ((cellhd->east > 180.0) && ((cellhd->east - 180.0) > epsilon_ew) 00117 && ((cellhd->east - 180.0) > GRASS_EPSILON)) { 00118 G_warning(_("Fixing subtle input data rounding error of east boundary (%g>%g)"), 00119 cellhd->east - 180.0, epsilon_ew); 00120 cellhd->east = 180.0; 00121 } 00122 #endif 00123 00124 while (cellhd->east <= cellhd->west) 00125 cellhd->east += 360.0; 00126 } 00127 00128 /* check the edge values */ 00129 if (cellhd->north <= cellhd->south) { 00130 if (cellhd->proj == PROJECTION_LL) 00131 return (_("North must be north of South")); 00132 else 00133 return (_("North must be larger than South")); 00134 } 00135 if (cellhd->east <= cellhd->west) 00136 return (_("East must be larger than West")); 00137 00138 /* compute rows and columns, if not set */ 00139 if (!row_flag) { 00140 cellhd->rows = 00141 (cellhd->north - cellhd->south + 00142 cellhd->ns_res / 2.0) / cellhd->ns_res; 00143 if (cellhd->rows == 0) 00144 cellhd->rows = 1; 00145 } 00146 if (!col_flag) { 00147 cellhd->cols = 00148 (cellhd->east - cellhd->west + 00149 cellhd->ew_res / 2.0) / cellhd->ew_res; 00150 if (cellhd->cols == 0) 00151 cellhd->cols = 1; 00152 } 00153 00154 if (cellhd->cols < 0 || cellhd->rows < 0) { 00155 return (_("Invalid coordinates")); 00156 } 00157 00158 00159 /* (re)compute the resolutions */ 00160 cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows; 00161 cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols; 00162 00163 return NULL; 00164 } 00165 00166 00194 char *G_adjust_Cell_head3(struct Cell_head *cellhd, int row_flag, 00195 int col_flag, int depth_flag) 00196 { 00197 if (!row_flag) { 00198 if (cellhd->ns_res <= 0) 00199 return (_("Illegal n-s resolution value")); 00200 if (cellhd->ns_res3 <= 0) 00201 return (_("Illegal n-s3 resolution value")); 00202 } 00203 else { 00204 if (cellhd->rows <= 0) 00205 return (_("Illegal row value")); 00206 if (cellhd->rows3 <= 0) 00207 return (_("Illegal row3 value")); 00208 } 00209 if (!col_flag) { 00210 if (cellhd->ew_res <= 0) 00211 return (_("Illegal e-w resolution value")); 00212 if (cellhd->ew_res3 <= 0) 00213 return (_("Illegal e-w3 resolution value")); 00214 } 00215 else { 00216 if (cellhd->cols <= 0) 00217 return (_("Illegal col value")); 00218 if (cellhd->cols3 <= 0) 00219 return (_("Illegal col3 value")); 00220 } 00221 if (!depth_flag) { 00222 if (cellhd->tb_res <= 0) 00223 return (_("Illegal t-b3 resolution value")); 00224 } 00225 else { 00226 if (cellhd->depths <= 0) 00227 return (_("Illegal depths value")); 00228 } 00229 00230 /* for lat/lon, check north,south. force east larger than west */ 00231 if (cellhd->proj == PROJECTION_LL) { 00232 double epsilon_ns, epsilon_ew; 00233 00234 /* TODO: find good thresholds */ 00235 epsilon_ns = 1. / cellhd->rows * 0.001; 00236 epsilon_ew = .000001; /* epsilon_ew calculation doesn't work due to cellhd->cols update/global wraparound below */ 00237 00238 G_debug(3, "G_adjust_Cell_head: epsilon_ns: %g, epsilon_ew: %g", 00239 epsilon_ns, epsilon_ew); 00240 00241 /* TODO: once working, change below G_warning to G_debug */ 00242 00243 /* fix rounding problems if input map slightly exceeds the world definition -180 90 180 -90 */ 00244 if (cellhd->north > 90.0) { 00245 if (((cellhd->north - 90.0) < epsilon_ns) && 00246 ((cellhd->north - 90.0) > GRASS_EPSILON)) { 00247 G_warning(_("Fixing subtle input data rounding error of north boundary (%g>%g)"), 00248 cellhd->north - 90.0, epsilon_ns); 00249 cellhd->north = 90.0; 00250 } 00251 else 00252 return (_("Illegal latitude for North")); 00253 } 00254 00255 if (cellhd->south < -90.0) { 00256 if (((cellhd->south + 90.0) < epsilon_ns) && 00257 ((cellhd->south + 90.0) < GRASS_EPSILON)) { 00258 G_warning(_("Fixing subtle input data rounding error of south boundary (%g>%g)"), 00259 cellhd->south + 90.0, epsilon_ns); 00260 cellhd->south = -90.0; 00261 } 00262 else 00263 return (_("Illegal latitude for South")); 00264 } 00265 00266 #if 0 00267 /* DISABLED: this breaks global wrap-around */ 00268 00269 G_debug(3, 00270 "G_adjust_Cell_head3() cellhd->west: %f, devi: %g, eps: %g", 00271 cellhd->west, cellhd->west + 180.0, epsilon_ew); 00272 00273 if ((cellhd->west < -180.0) && ((cellhd->west + 180.0) < epsilon_ew) 00274 && ((cellhd->west + 180.0) < GRASS_EPSILON)) { 00275 G_warning(_("Fixing subtle input data rounding error of west boundary (%g>%g)"), 00276 cellhd->west + 180.0, epsilon_ew); 00277 cellhd->west = -180.0; 00278 } 00279 00280 G_debug(3, 00281 "G_adjust_Cell_head3() cellhd->east: %f, devi: %g, eps: %g", 00282 cellhd->east, cellhd->east - 180.0, epsilon_ew); 00283 00284 if ((cellhd->east > 180.0) && ((cellhd->east - 180.0) > epsilon_ew) 00285 && ((cellhd->east - 180.0) > GRASS_EPSILON)) { 00286 G_warning(_("Fixing subtle input data rounding error of east boundary (%g>%g)"), 00287 cellhd->east - 180.0, epsilon_ew); 00288 cellhd->east = 180.0; 00289 } 00290 #endif 00291 00292 while (cellhd->east <= cellhd->west) 00293 cellhd->east += 360.0; 00294 } 00295 00296 /* check the edge values */ 00297 if (cellhd->north <= cellhd->south) { 00298 if (cellhd->proj == PROJECTION_LL) 00299 return (_("North must be north of South")); 00300 else 00301 return (_("North must be larger than South")); 00302 } 00303 if (cellhd->east <= cellhd->west) 00304 return (_("East must be larger than West")); 00305 if (cellhd->top <= cellhd->bottom) 00306 return (_("Top must be larger than Bottom")); 00307 00308 00309 /* compute rows and columns, if not set */ 00310 if (!row_flag) { 00311 cellhd->rows = 00312 (cellhd->north - cellhd->south + 00313 cellhd->ns_res / 2.0) / cellhd->ns_res; 00314 if (cellhd->rows == 0) 00315 cellhd->rows = 1; 00316 00317 cellhd->rows3 = 00318 (cellhd->north - cellhd->south + 00319 cellhd->ns_res3 / 2.0) / cellhd->ns_res3; 00320 if (cellhd->rows3 == 0) 00321 cellhd->rows3 = 1; 00322 } 00323 if (!col_flag) { 00324 cellhd->cols = 00325 (cellhd->east - cellhd->west + 00326 cellhd->ew_res / 2.0) / cellhd->ew_res; 00327 if (cellhd->cols == 0) 00328 cellhd->cols = 1; 00329 00330 cellhd->cols3 = 00331 (cellhd->east - cellhd->west + 00332 cellhd->ew_res3 / 2.0) / cellhd->ew_res3; 00333 if (cellhd->cols3 == 0) 00334 cellhd->cols3 = 1; 00335 } 00336 00337 if (!depth_flag) { 00338 cellhd->depths = 00339 (cellhd->top - cellhd->bottom + 00340 cellhd->tb_res / 2.0) / cellhd->tb_res; 00341 if (cellhd->depths == 0) 00342 cellhd->depths = 1; 00343 00344 } 00345 00346 if (cellhd->cols < 0 || cellhd->rows < 0 || cellhd->cols3 < 0 || 00347 cellhd->rows3 < 0 || cellhd->depths < 0) { 00348 return (_("Invalid coordinates")); 00349 } 00350 00351 /* (re)compute the resolutions */ 00352 cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows; 00353 cellhd->ns_res3 = (cellhd->north - cellhd->south) / cellhd->rows3; 00354 cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols; 00355 cellhd->ew_res3 = (cellhd->east - cellhd->west) / cellhd->cols3; 00356 cellhd->tb_res = (cellhd->top - cellhd->bottom) / cellhd->depths; 00357 00358 return NULL; 00359 }