adj_cellhd.c

Go to the documentation of this file.
00001 #include <grass/gis.h>
00002 #include <grass/glocale.h>
00003 
00004 /***********************************************************
00005 
00006 char *
00007 G_adjust_Cell_head (cellhd, row_flag, col_flag)
00008     struct Cell_head *cellhd;
00009 
00010  This function fills in missing parts of the input cell header.
00011  It also makes projection-specific adjustments (such as lat-lon).
00012 
00013  The north, south, east, west parts of the header must be set
00014  before calling this routine.
00015 
00016  The other parts depend on row_flag and col_flag:
00017 
00018  row_flag TRUE
00019     the rows in the cell header is already set
00020     the n-s resolution is to be computed.
00021 
00022  row_flag FALSE
00023     the n-s resolution in the cell header is already set,
00024     the rows are to be computed.
00025     (the resolution may be adjusted to conform to the computed rows)
00026 
00027  col_flag TRUE
00028     the cols in the cell header is already set
00029     the e-w resolution is to be computed.
00030 
00031  col_flag FALSE
00032     the e-w resolution in the cell header is already set,
00033     the cols are to be computed.
00034     (the resolution may be adjusted to conform to the computed cols)
00035 
00036  returns an error message, or NULL if ok
00037  note: don't free this error message.
00038 ************************************************************/
00039 
00064 char *G_adjust_Cell_head(struct Cell_head *cellhd,int row_flag,int col_flag)
00065 {
00066     if (!row_flag)
00067     {
00068         if (cellhd->ns_res <= 0)
00069             return (_("Illegal n-s resolution value"));
00070     }
00071     else
00072     {
00073         if (cellhd->rows <= 0)
00074             return (_("Illegal row value"));
00075     }
00076     if (!col_flag)
00077     {
00078         if (cellhd->ew_res <= 0)
00079             return (_("Illegal e-w resolution value"));
00080     }
00081     else
00082     {
00083         if (cellhd->cols <= 0)
00084             return (_("Illegal col value"));
00085     }
00086 
00087 /* for lat/lon, check north,south. force east larger than west */
00088     if (cellhd->proj == PROJECTION_LL)
00089     {
00090         double epsilon_ns, epsilon_ew;
00091         
00092         /* TODO: find good thresholds */
00093         epsilon_ns= 1. / cellhd->rows * 0.001;
00094         epsilon_ew= .000001;  /* epsilon_ew calculation doesn't work due to cellhd->cols update/global wraparound below */
00095 
00096         G_debug(3,"G_adjust_Cell_head: epsilon_ns: %g, epsilon_ew: %g", epsilon_ns, epsilon_ew);
00097         
00098         /* TODO: once working, change below G_warning to G_debug */
00099         
00100         /* fix rounding problems if input map slightly exceeds the world definition -180 90 180 -90 */
00101         if (cellhd->north > 90.0) {
00102             if ( ((cellhd->north - 90.0) < epsilon_ns) && ((cellhd->north - 90.0) > GRASS_EPSILON) ) {
00103                 G_warning (_("Fixing subtle input data rounding error of north boundary (%g>%g)"), cellhd->north - 90.0, epsilon_ns);
00104                 cellhd->north = 90.0;
00105             } else
00106                 return (_("Illegal latitude for North"));
00107         }
00108 
00109         if (cellhd->south < -90.0) {
00110             if ( ((cellhd->south + 90.0) < epsilon_ns) && ((cellhd->south + 90.0) < GRASS_EPSILON) ) {
00111                 G_warning (_("Fixing subtle input data rounding error of south boundary (%g>%g)"), cellhd->south + 90.0, epsilon_ns);
00112                 cellhd->south = -90.0;
00113             } else
00114                 return (_("Illegal latitude for South"));
00115         }
00116 
00117 /* DISABLED: this breaks global wrap-around
00118         G_debug(3,"G_adjust_Cell_head()  cellhd->west: %f, devi: %g, eps: %g", 
00119            cellhd->west, cellhd->west + 180.0, epsilon_ew);
00120         if ( (cellhd->west < -180.0) && ((cellhd->west + 180.0) < epsilon_ew)
00121            && ((cellhd->west + 180.0) < GRASS_EPSILON) ) {
00122             G_warning (_("Fixing subtle input data rounding error of west boundary (%g>%g)"),
00123                cellhd->west + 180.0, epsilon_ew);
00124             cellhd->west = -180.0;
00125         }
00126 
00127         G_debug(3,"G_adjust_Cell_head()  cellhd->east: %f, devi: %g, eps: %g",
00128            cellhd->east, cellhd->east - 180.0, epsilon_ew);
00129         if ( (cellhd->east > 180.0) && ((cellhd->east - 180.0) > epsilon_ew)
00130            && ((cellhd->east - 180.0) > GRASS_EPSILON) ) {
00131             G_warning (_("Fixing subtle input data rounding error of east boundary (%g>%g)"),
00132                cellhd->east - 180.0, epsilon_ew);
00133             cellhd->east = 180.0;
00134         }
00135 */
00136         while (cellhd->east <= cellhd->west)
00137             cellhd->east += 360.0;
00138     }
00139 
00140 /* check the edge values */
00141     if (cellhd->north <= cellhd->south)
00142     {
00143         if (cellhd->proj == PROJECTION_LL)
00144             return (_("North must be north of South"));
00145         else
00146             return (_("North must be larger than South"));
00147     }
00148     if (cellhd->east <= cellhd->west)
00149         return (_("East must be larger than West"));
00150 
00151 
00152 /* compute rows and columns, if not set */
00153     if (!row_flag)
00154     {
00155         cellhd->rows = (cellhd->north - cellhd->south + cellhd->ns_res/2.0) / cellhd->ns_res;
00156         if (cellhd->rows == 0)
00157             cellhd->rows = 1;
00158     }
00159     if (!col_flag)
00160     {
00161         cellhd->cols = (cellhd->east - cellhd->west + cellhd->ew_res/2.0) / cellhd->ew_res;
00162         if (cellhd->cols == 0)
00163             cellhd->cols = 1;
00164     }
00165 
00166     if (cellhd->cols < 0 || cellhd->rows < 0 )
00167     {
00168         return (_("Invalid coordinates"));
00169     }
00170 
00171 
00172 /* (re)compute the resolutions */
00173     cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
00174     cellhd->ew_res = (cellhd->east  - cellhd->west)  / cellhd->cols;
00175         
00176     return NULL;
00177 }
00178 
00203 char *G_adjust_Cell_head3(struct Cell_head *cellhd,int row_flag,int col_flag, int depth_flag)
00204 {
00205     if (!row_flag)
00206     {
00207         if (cellhd->ns_res <= 0)
00208             return (_("Illegal n-s resolution value"));
00209         if (cellhd->ns_res3 <= 0)
00210             return (_("Illegal n-s3 resolution value"));
00211     }
00212     else
00213     {
00214         if (cellhd->rows <= 0)
00215             return (_("Illegal row value"));
00216         if (cellhd->rows3 <= 0)
00217             return (_("Illegal row3 value"));
00218     }
00219     if (!col_flag)
00220     {
00221         if (cellhd->ew_res <= 0)
00222             return (_("Illegal e-w resolution value"));
00223         if (cellhd->ew_res3 <= 0)
00224             return (_("Illegal e-w3 resolution value"));
00225     }
00226     else
00227     {
00228         if (cellhd->cols <= 0)
00229             return (_("Illegal col value"));
00230         if (cellhd->cols3 <= 0)
00231             return (_("Illegal col3 value"));
00232     }
00233     if (!depth_flag)
00234     {
00235         if (cellhd->tb_res <= 0)
00236             return (_("Illegal t-b3 resolution value"));
00237     }
00238     else
00239     {
00240         if (cellhd->depths <= 0)
00241             return (_("Illegal depths value"));
00242     }
00243 
00244 /* for lat/lon, check north,south. force east larger than west */
00245     if (cellhd->proj == PROJECTION_LL)
00246     {
00247         double epsilon_ns, epsilon_ew;
00248         
00249         /* TODO: find good thresholds */
00250         epsilon_ns= 1. / cellhd->rows * 0.001;
00251         epsilon_ew= .000001;  /* epsilon_ew calculation doesn't work due to cellhd->cols update/global wraparound below */
00252 
00253         G_debug(3,"G_adjust_Cell_head: epsilon_ns: %g, epsilon_ew: %g", epsilon_ns, epsilon_ew);
00254         
00255         /* TODO: once working, change below G_warning to G_debug */
00256         
00257         /* fix rounding problems if input map slightly exceeds the world definition -180 90 180 -90 */
00258         if (cellhd->north > 90.0) {
00259             if ( ((cellhd->north - 90.0) < epsilon_ns) && ((cellhd->north - 90.0) > GRASS_EPSILON) ) {
00260                 G_warning (_("Fixing subtle input data rounding error of north boundary (%g>%g)"), cellhd->north - 90.0, epsilon_ns);
00261                 cellhd->north = 90.0;
00262             } else
00263                 return (_("Illegal latitude for North"));
00264         }
00265 
00266         if (cellhd->south < -90.0) {
00267             if ( ((cellhd->south + 90.0) < epsilon_ns) && ((cellhd->south + 90.0) < GRASS_EPSILON) ) {
00268                 G_warning (_("Fixing subtle input data rounding error of south boundary (%g>%g)"), cellhd->south + 90.0, epsilon_ns);
00269                 cellhd->south = -90.0;
00270             } else
00271                 return (_("Illegal latitude for South"));
00272         }
00273 
00274 /* DISABLED: this breaks global wrap-around
00275         G_debug(3,"G_adjust_Cell_head3() cellhd->west: %f, devi: %g, eps: %g",
00276            cellhd->west, cellhd->west + 180.0, epsilon_ew);
00277         if ( (cellhd->west < -180.0) && ((cellhd->west + 180.0) < epsilon_ew)
00278            && ((cellhd->west + 180.0) < GRASS_EPSILON) ) {
00279             G_warning (_("Fixing subtle input data rounding error of west boundary (%g>%g)"),
00280                cellhd->west + 180.0, epsilon_ew);
00281             cellhd->west = -180.0;
00282         }
00283 
00284         G_debug(3,"G_adjust_Cell_head3() cellhd->east: %f, devi: %g, eps: %g",
00285            cellhd->east, cellhd->east - 180.0, epsilon_ew);
00286         if ( (cellhd->east > 180.0) && ((cellhd->east - 180.0) > epsilon_ew)
00287            && ((cellhd->east - 180.0) > GRASS_EPSILON) ) {
00288             G_warning (_("Fixing subtle input data rounding error of east boundary (%g>%g)"),
00289                cellhd->east - 180.0, epsilon_ew);
00290             cellhd->east = 180.0;
00291         }
00292 */
00293         while (cellhd->east <= cellhd->west)
00294             cellhd->east += 360.0;
00295     }
00296 
00297 /* check the edge values */
00298     if (cellhd->north <= cellhd->south)
00299     {
00300         if (cellhd->proj == PROJECTION_LL)
00301             return (_("North must be north of South"));
00302         else
00303             return (_("North must be larger than South"));
00304     }
00305     if (cellhd->east <= cellhd->west)
00306         return (_("East must be larger than West"));
00307     if (cellhd->top <= cellhd->bottom)
00308         return (_("Top must be larger than Bottom"));
00309     
00310 
00311 /* compute rows and columns, if not set */
00312     if (!row_flag)
00313     {
00314         cellhd->rows = (cellhd->north - cellhd->south + cellhd->ns_res/2.0) / cellhd->ns_res;
00315         if (cellhd->rows == 0)
00316             cellhd->rows = 1;
00317 
00318         cellhd->rows3 = (cellhd->north - cellhd->south + cellhd->ns_res3/2.0) / cellhd->ns_res3;
00319         if (cellhd->rows3 == 0)
00320             cellhd->rows3 = 1;
00321     }
00322     if (!col_flag)
00323     {
00324         cellhd->cols = (cellhd->east - cellhd->west + cellhd->ew_res/2.0) / cellhd->ew_res;
00325         if (cellhd->cols == 0)
00326             cellhd->cols = 1;
00327 
00328         cellhd->cols3 = (cellhd->east - cellhd->west + cellhd->ew_res3/2.0) / cellhd->ew_res3;
00329         if (cellhd->cols3 == 0)
00330             cellhd->cols3 = 1;
00331     }
00332 
00333     if (!depth_flag)
00334     {
00335         cellhd->depths = (cellhd->top - cellhd->bottom + cellhd->tb_res/2.0) / cellhd->tb_res;
00336         if (cellhd->depths == 0)
00337             cellhd->depths = 1;
00338 
00339     }
00340 
00341 
00342     if (cellhd->cols < 0 || cellhd->rows < 0 || cellhd->cols3 < 0 || cellhd->rows3 < 0 
00343         || cellhd->depths < 0 )
00344     {
00345         return (_("Invalid coordinates"));
00346     }
00347 
00348 
00349 /* (re)compute the resolutions */
00350     cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
00351     cellhd->ns_res3 = (cellhd->north - cellhd->south) / cellhd->rows3;
00352     cellhd->ew_res = (cellhd->east  - cellhd->west)  / cellhd->cols;
00353     cellhd->ew_res3 = (cellhd->east  - cellhd->west)  / cellhd->cols3;
00354     cellhd->tb_res = (cellhd->top  - cellhd->bottom)  / cellhd->depths;
00355         
00356     return NULL;
00357 }
00358 

Generated on Wed Dec 19 14:59:05 2007 for GRASS by  doxygen 1.5.4