GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71769
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
adj_cellhd.c
Go to the documentation of this file.
1 /*!
2  * \file lib/gis/adj_cellhd.c
3  *
4  * \brief GIS Library - CELL header adjustment.
5  *
6  * (C) 2001-2009 by the GRASS Development Team
7  *
8  * This program is free software under the GNU General Public License
9  * (>=v2). Read the file COPYING that comes with GRASS for details.
10  *
11  * \author Original author CERL
12  */
13 
14 #include <math.h>
15 #include <string.h>
16 #include <grass/gis.h>
17 #include <grass/glocale.h>
18 
19 #define LL_TOLERANCE 10
20 
21 /* TODO: find good thresholds */
22 /* deviation measured in cells */
23 static double llepsilon = 0.01;
24 static double fpepsilon = 1.0e-9;
25 
26 static int ll_wrap(struct Cell_head *cellhd);
27 static int ll_check_ns(struct Cell_head *cellhd);
28 static int ll_check_ew(struct Cell_head *cellhd);
29 
30 /*!
31  * \brief Adjust cell header.
32  *
33  * This function fills in missing parts of the input cell header (or
34  * region). It also makes projection-specific adjustments. The
35  * <i>cellhd</i> structure must have its <i>north, south, east,
36  * west</i>, and <i>proj</i> fields set.
37  *
38  * If <i>row_flag</i> is true, then the north-south resolution is
39  * computed from the number of <i>rows</i> in the <i>cellhd</i>
40  * structure. Otherwise the number of <i>rows</i> is computed from the
41  * north-south resolution in the structure, similarly for
42  * <i>col_flag</i> and the number of columns and the east-west
43  * resolution.
44  *
45  * <b>Note:</b> 3D values are not adjusted.
46  *
47  * \param[in,out] cellhd pointer to Cell_head structure
48  * \param row_flag compute n-s resolution
49  * \param col_flag compute e-w resolution
50  */
51 void G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag)
52 {
53  double old_res;
54 
55  if (!row_flag) {
56  if (cellhd->ns_res <= 0)
57  G_fatal_error(_("Illegal n-s resolution value <%lf>"), cellhd->ns_res);
58  }
59  else {
60  if (cellhd->rows <= 0)
61  G_fatal_error(_("Illegal row value"));
62  }
63  if (!col_flag) {
64  if (cellhd->ew_res <= 0)
65  G_fatal_error(_("Illegal e-w resolution value"));
66  }
67  else {
68  if (cellhd->cols <= 0)
69  G_fatal_error(_("Illegal col value"));
70  }
71 
72  /* check the edge values */
73  if (cellhd->north <= cellhd->south) {
74  if (cellhd->proj == PROJECTION_LL)
75  G_fatal_error(_("North must be north of South"));
76  else
77  G_fatal_error(_("North must be larger than South"));
78  }
79 
80  ll_wrap(cellhd);
81 
82  if (cellhd->east <= cellhd->west)
83  G_fatal_error(_("East must be larger than West"));
84 
85  /* compute rows and columns, if not set */
86  if (!row_flag) {
87  cellhd->rows =
88  (cellhd->north - cellhd->south +
89  cellhd->ns_res / 2.0) / cellhd->ns_res;
90  if (cellhd->rows == 0)
91  cellhd->rows = 1;
92  }
93  if (!col_flag) {
94  cellhd->cols =
95  (cellhd->east - cellhd->west +
96  cellhd->ew_res / 2.0) / cellhd->ew_res;
97  if (cellhd->cols == 0)
98  cellhd->cols = 1;
99  }
100 
101  if (cellhd->cols < 0 || cellhd->rows < 0) {
102  G_fatal_error(_("Invalid coordinates"));
103  }
104 
105  /* (re)compute the resolutions */
106  old_res = cellhd->ns_res;
107  cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
108  if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
109  G_verbose_message(_("NS resolution has been changed"));
110 
111  old_res = cellhd->ew_res;
112  cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
113  if (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
114  G_verbose_message(_("EW resolution has been changed"));
115 
116  if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
117  G_verbose_message(_("NS and EW resolutions are different"));
118 
119  ll_check_ns(cellhd);
120  ll_check_ew(cellhd);
121 }
122 
123 /*!
124  * \brief Adjust cell header for 3D values.
125  *
126  * This function fills in missing parts of the input cell header (or
127  * region). It also makes projection-specific adjustments. The
128  * <i>cellhd</i> structure must have its <i>north, south, east,
129  * west</i>, and <i>proj</i> fields set.
130  *
131  * If <i>row_flag</i> is true, then the north-south resolution is computed
132  * from the number of <i>rows</i> in the <i>cellhd</i> structure.
133  * Otherwise the number of <i>rows</i> is computed from the north-south
134  * resolution in the structure, similarly for <i>col_flag</i> and the
135  * number of columns and the east-west resolution.
136  *
137  * If <i>depth_flag</i> is true, top-bottom resolution is calculated
138  * from depths.
139  * If <i>depth_flag</i> are false, number of depths is calculated from
140  * top-bottom resolution.
141  *
142  * \warning This function can cause segmentation fault without any warning
143  * when it is called with Cell_head top and bottom set to zero.
144  *
145  * \param[in,out] cellhd pointer to Cell_head structure
146  * \param row_flag compute n-s resolution
147  * \param col_flag compute e-w resolution
148  * \param depth_flag compute t-b resolution
149  */
150 void G_adjust_Cell_head3(struct Cell_head *cellhd, int row_flag,
151  int col_flag, int depth_flag)
152 {
153  double old_res;
154 
155  if (!row_flag) {
156  if (cellhd->ns_res <= 0)
157  G_fatal_error(_("Illegal n-s resolution value"));
158  if (cellhd->ns_res3 <= 0)
159  G_fatal_error(_("Illegal n-s3 resolution value"));
160  }
161  else {
162  if (cellhd->rows <= 0)
163  G_fatal_error(_("Illegal row value"));
164  if (cellhd->rows3 <= 0)
165  G_fatal_error(_("Illegal row3 value"));
166  }
167  if (!col_flag) {
168  if (cellhd->ew_res <= 0)
169  G_fatal_error(_("Illegal e-w resolution value"));
170  if (cellhd->ew_res3 <= 0)
171  G_fatal_error(_("Illegal e-w3 resolution value"));
172  }
173  else {
174  if (cellhd->cols <= 0)
175  G_fatal_error(_("Illegal col value"));
176  if (cellhd->cols3 <= 0)
177  G_fatal_error(_("Illegal col3 value"));
178  }
179  if (!depth_flag) {
180  if (cellhd->tb_res <= 0)
181  G_fatal_error(_("Illegal t-b3 resolution value"));
182  }
183  else {
184  if (cellhd->depths <= 0)
185  G_fatal_error(_("Illegal depths value"));
186  }
187 
188  /* check the edge values */
189  if (cellhd->north <= cellhd->south) {
190  if (cellhd->proj == PROJECTION_LL)
191  G_fatal_error(_("North must be north of South"));
192  else
193  G_fatal_error(_("North must be larger than South"));
194  }
195 
196  ll_wrap(cellhd);
197 
198  if (cellhd->east <= cellhd->west)
199  G_fatal_error(_("East must be larger than West"));
200 
201  if (cellhd->top <= cellhd->bottom)
202  G_fatal_error(_("Top must be larger than Bottom"));
203 
204  /* compute rows and columns, if not set */
205  if (!row_flag) {
206  cellhd->rows =
207  (cellhd->north - cellhd->south +
208  cellhd->ns_res / 2.0) / cellhd->ns_res;
209  if (cellhd->rows == 0)
210  cellhd->rows = 1;
211 
212  cellhd->rows3 =
213  (cellhd->north - cellhd->south +
214  cellhd->ns_res3 / 2.0) / cellhd->ns_res3;
215  if (cellhd->rows3 == 0)
216  cellhd->rows3 = 1;
217  }
218  if (!col_flag) {
219  cellhd->cols =
220  (cellhd->east - cellhd->west +
221  cellhd->ew_res / 2.0) / cellhd->ew_res;
222  if (cellhd->cols == 0)
223  cellhd->cols = 1;
224 
225  cellhd->cols3 =
226  (cellhd->east - cellhd->west +
227  cellhd->ew_res3 / 2.0) / cellhd->ew_res3;
228  if (cellhd->cols3 == 0)
229  cellhd->cols3 = 1;
230  }
231 
232  if (!depth_flag) {
233  cellhd->depths =
234  (cellhd->top - cellhd->bottom +
235  cellhd->tb_res / 2.0) / cellhd->tb_res;
236  if (cellhd->depths == 0)
237  cellhd->depths = 1;
238  }
239 
240  if (cellhd->cols < 0 || cellhd->rows < 0 || cellhd->cols3 < 0 ||
241  cellhd->rows3 < 0 || cellhd->depths < 0) {
242  G_fatal_error(_("Invalid coordinates"));
243  }
244 
245  /* (re)compute the resolutions */
246  old_res = cellhd->ns_res;
247  cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
248  if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
249  G_verbose_message(_("NS resolution has been changed"));
250 
251  old_res = cellhd->ew_res;
252  cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
253  if (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
254  G_verbose_message(_("EW resolution has been changed"));
255 
256  if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
257  G_verbose_message(_("NS and EW resolutions are different"));
258 
259  ll_check_ns(cellhd);
260  ll_check_ew(cellhd);
261 
262  cellhd->ns_res3 = (cellhd->north - cellhd->south) / cellhd->rows3;
263  cellhd->ew_res3 = (cellhd->east - cellhd->west) / cellhd->cols3;
264  cellhd->tb_res = (cellhd->top - cellhd->bottom) / cellhd->depths;
265 }
266 
267 static int ll_wrap(struct Cell_head *cellhd)
268 {
269  double shift;
270 
271  /* for lat/lon, force east larger than west, try to wrap to -180, 180 */
272  if (cellhd->proj != PROJECTION_LL)
273  return 0;
274 
275  if (cellhd->east <= cellhd->west) {
276  G_warning(_("East (%.15g) is not larger than West (%.15g)"),
277  cellhd->east, cellhd->west);
278 
279  while (cellhd->east <= cellhd->west)
280  cellhd->east += 360.0;
281  }
282 
283  /* with east larger than west,
284  * any 360 degree W-E extent can be represented within -360, 360
285  * but not within -180, 180 */
286 
287  /* try to shift to within -180, 180 */
288  shift = 0;
289  while (cellhd->west + shift >= 180) {
290  shift -= 360.0;
291  }
292  while (cellhd->east + shift <= -180) {
293  shift += 360.0;
294  }
295 
296  /* try to shift to within -360, 360 */
297  while (cellhd->east + shift > 360) {
298  shift -= 360.0;
299  }
300  while (cellhd->west + shift <= -360) {
301  shift += 360.0;
302  }
303 
304  if (shift) {
305  cellhd->west += shift;
306  cellhd->east += shift;
307  }
308 
309  /* very liberal thresholds */
310  if (cellhd->north > 90.0 + LL_TOLERANCE)
311  G_fatal_error(_("Illegal latitude for North: %g"), cellhd->north);
312  if (cellhd->south < -90.0 - LL_TOLERANCE)
313  G_fatal_error(_("Illegal latitude for South: %g"), cellhd->south);
314 
315 #if 0
316  /* disabled: allow W-E extents larger than 360 degree e.g. for display */
317  if (cellhd->west < -360.0 - LL_TOLERANCE) {
318  G_debug(1, "East: %g", cellhd->east);
319  G_fatal_error(_("Illegal longitude for West: %g"), cellhd->west);
320  }
321  if (cellhd->east > 360.0 + LL_TOLERANCE) {
322  G_debug(1, "West: %g", cellhd->west);
323  G_fatal_error(_("Illegal longitude for East: %g"), cellhd->east);
324  }
325 #endif
326 
327  return 1;
328 }
329 
330 static int ll_check_ns(struct Cell_head *cellhd)
331 {
332  int lladjust;
333  double diff;
334  int ncells;
335 
336  /* lat/lon checks */
337  if (cellhd->proj != PROJECTION_LL)
338  return 0;
339 
340  lladjust = 0;
341 
342  G_debug(3, "ll_check_ns: epsilon: %g", llepsilon);
343 
344  /* North, South: allow a half cell spill-over */
345 
346  diff = (cellhd->north - cellhd->south) / cellhd->ns_res;
347  ncells = (int) (diff + 0.5);
348  diff -= ncells;
349  if ((diff < 0 && diff < -fpepsilon) ||
350  (diff > 0 && diff > fpepsilon)) {
351  G_verbose_message(_("NS extent does not match NS resolution: %g cells difference"),
352  diff);
353  }
354 
355  /* north */
356  diff = (cellhd->north - 90) / cellhd->ns_res;
357  if (diff < 0)
358  diff = -diff;
359  if (cellhd->north < 90.0 && diff < 1.0 ) {
360  G_verbose_message(_("%g cells missing to reach 90 degree north"),
361  diff);
362  if (diff < llepsilon && diff > fpepsilon) {
363  G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
364  cellhd->north - 90.0);
365  /* check only, do not modify
366  cellhd->north = 90.0;
367  lladjust = 1;
368  */
369  }
370  }
371  if (cellhd->north > 90.0) {
372  if (diff <= 0.5 + llepsilon) {
373  G_important_message(_("90 degree north is exceeded by %g cells"),
374  diff);
375 
376  if (diff < llepsilon && diff > fpepsilon) {
377  G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
378  cellhd->north - 90.0);
379  G_debug(1, "North of north in seconds: %g",
380  (cellhd->north - 90.0) * 3600);
381  /* check only, do not modify
382  cellhd->north = 90.0;
383  lladjust = 1;
384  */
385  }
386 
387  diff = diff - 0.5;
388  if (diff < 0)
389  diff = -diff;
390  if (diff < llepsilon && diff > fpepsilon) {
391  G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
392  cellhd->north - 90.0 - cellhd->ns_res / 2.0);
393  G_debug(1, "North of north + 0.5 cells in seconds: %g",
394  (cellhd->north - 90.0 - cellhd->ns_res / 2.0) * 3600);
395  /* check only, do not modify
396  cellhd->north = 90.0 + cellhd->ns_res / 2.0;
397  lladjust = 1;
398  */
399  }
400  }
401  else
402  G_fatal_error(_("Illegal latitude for North"));
403  }
404 
405  /* south */
406  diff = (cellhd->south + 90) / cellhd->ns_res;
407  if (diff < 0)
408  diff = -diff;
409  if (cellhd->south > -90.0 && diff < 1.0 ) {
410  G_verbose_message(_("%g cells missing to reach 90 degree south"),
411  diff);
412  if (diff < llepsilon && diff > fpepsilon) {
413  G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
414  cellhd->south + 90.0);
415  /* check only, do not modify
416  cellhd->south = -90.0;
417  lladjust = 1;
418  */
419  }
420  }
421  if (cellhd->south < -90.0) {
422  if (diff <= 0.5 + llepsilon) {
423  G_important_message(_("90 degree south is exceeded by %g cells"),
424  diff);
425 
426  if (diff < llepsilon && diff > fpepsilon) {
427  G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
428  cellhd->south + 90);
429  G_debug(1, "South of south in seconds: %g",
430  (-cellhd->south - 90) * 3600);
431  /* check only, do not modify
432  cellhd->south = -90.0;
433  lladjust = 1;
434  */
435  }
436 
437  diff = diff - 0.5;
438  if (diff < 0)
439  diff = -diff;
440  if (diff < llepsilon && diff > fpepsilon) {
441  G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
442  cellhd->south + 90 + cellhd->ns_res / 2.0);
443  G_debug(1, "South of south + 0.5 cells in seconds: %g",
444  (-cellhd->south - 90 - cellhd->ns_res / 2.0) * 3600);
445  /* check only, do not modify
446  cellhd->south = -90.0 - cellhd->ns_res / 2.0;
447  lladjust = 1;
448  */
449  }
450  }
451  else
452  G_fatal_error(_("Illegal latitude for South"));
453  }
454 
455  if (lladjust)
456  cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
457 
458  return lladjust;
459 }
460 
461 static int ll_check_ew(struct Cell_head *cellhd)
462 {
463  int lladjust;
464  double diff;
465  int ncells;
466 
467  /* lat/lon checks */
468  if (cellhd->proj != PROJECTION_LL)
469  return 0;
470 
471  lladjust = 0;
472 
473  G_debug(3, "ll_check_ew: epsilon: %g", llepsilon);
474 
475  /* west - east, no adjustment */
476  diff = (cellhd->east - cellhd->west) / cellhd->ew_res;
477  ncells = (int) (diff + 0.5);
478  diff -= ncells;
479  if ((diff < 0 && diff < -fpepsilon) ||
480  (diff > 0 && diff > fpepsilon)) {
481  G_verbose_message(_("EW extent does not match EW resolution: %g cells difference"),
482  diff);
483  }
484  if (cellhd->east - cellhd->west > 360.0) {
485  diff = (cellhd->east - cellhd->west - 360.0) / cellhd->ew_res;
486  if (diff > fpepsilon)
487  G_important_message(_("360 degree EW extent is exceeded by %g cells"),
488  diff);
489  }
490  else if (cellhd->east - cellhd->west < 360.0) {
491  diff = (360.0 - (cellhd->east - cellhd->west)) / cellhd->ew_res;
492  if (diff < 1.0 && diff > fpepsilon)
493  G_verbose_message(_("%g cells missing to cover 360 degree EW extent"),
494  diff);
495  }
496 
497  return lladjust;
498 }
499 
500 /*!
501  * \brief Adjust window for lat/lon.
502  *
503  * This function tries to automatically fix fp precision issues and
504  * adjust rounding errors for lat/lon.
505  *
506  * <b>Note:</b> 3D values are not adjusted.
507  *
508  * \param[in,out] cellhd pointer to Cell_head structure
509  * \return 1 if window was adjusted
510  * \return 0 if window was not adjusted
511  */
512 int G_adjust_window_ll(struct Cell_head *cellhd)
513 {
514  int ll_adjust, res_adj;
515  double dsec, dsec2;
516  char buf[100], buf2[100];
517  double diff;
518  double old, new;
519  struct Cell_head cellhds; /* everyting in seconds, not degrees */
520 
521  /* lat/lon checks */
522  if (cellhd->proj != PROJECTION_LL)
523  return 0;
524 
525  /* put everything through ll_format + ll_scan */
526  G_llres_format(cellhd->ns_res, buf);
527  if (G_llres_scan(buf, &new) != 1)
528  G_fatal_error(_("Invalid NS resolution"));
529  cellhd->ns_res = new;
530 
531  G_llres_format(cellhd->ew_res, buf);
532  if (G_llres_scan(buf, &new) != 1)
533  G_fatal_error(_("Invalid EW resolution"));
534  cellhd->ew_res = new;
535 
536  G_lat_format(cellhd->north, buf);
537  if (G_lat_scan(buf, &new) != 1)
538  G_fatal_error(_("Invalid North"));
539  cellhd->north = new;
540 
541  G_lat_format(cellhd->south, buf);
542  if (G_lat_scan(buf, &new) != 1)
543  G_fatal_error(_("Invalid South"));
544  cellhd->south = new;
545 
546  G_lon_format(cellhd->west, buf);
547  if (G_lon_scan(buf, &new) != 1)
548  G_fatal_error(_("Invalid West"));
549  cellhd->west = new;
550 
551  G_lon_format(cellhd->east, buf);
552  if (G_lon_scan(buf, &new) != 1)
553  G_fatal_error(_("Invalid East"));
554  cellhd->east = new;
555 
556  /* convert to seconds */
557  cellhds = *cellhd;
558 
559  old = cellhds.ns_res * 3600;
560  sprintf(buf, "%f", old);
561  sscanf(buf, "%lf", &new);
562  cellhds.ns_res = new;
563 
564  old = cellhds.ew_res * 3600;
565  sprintf(buf, "%f", old);
566  sscanf(buf, "%lf", &new);
567  cellhds.ew_res = new;
568 
569  old = cellhds.north * 3600;
570  sprintf(buf, "%f", old);
571  sscanf(buf, "%lf", &new);
572  cellhds.north = new;
573 
574  old = cellhds.south * 3600;
575  sprintf(buf, "%f", old);
576  sscanf(buf, "%lf", &new);
577  cellhds.south = new;
578 
579  old = cellhds.west * 3600;
580  sprintf(buf, "%f", old);
581  sscanf(buf, "%lf", &new);
582  cellhds.west = new;
583 
584  old = cellhds.east * 3600;
585  sprintf(buf, "%f", old);
586  sscanf(buf, "%lf", &new);
587  cellhds.east = new;
588 
589  ll_adjust = 0;
590 
591  /* N - S */
592  /* resolution */
593  res_adj = 0;
594  old = cellhds.ns_res;
595 
596  if (old > 0.4) {
597  /* round to nearest 0.1 sec */
598  dsec = old * 10;
599  dsec2 = floor(dsec + 0.5);
600  new = dsec2 / 10;
601  diff = fabs(dsec2 - dsec) / dsec;
602  if (diff > 0 && diff < llepsilon) {
603  G_llres_format(old / 3600, buf);
604  G_llres_format(new / 3600, buf2);
605  if (strcmp(buf, buf2))
606  G_verbose_message(_("NS resolution rounded from %s to %s"),
607  buf, buf2);
608  ll_adjust = 1;
609  res_adj = 1;
610  cellhds.ns_res = new;
611  }
612  }
613 
614  if (res_adj) {
615  double n_off, s_off;
616 
617  old = cellhds.north;
618  dsec = old * 10;
619  dsec2 = floor(dsec + 0.5);
620  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
621  n_off = diff;
622 
623  old = cellhds.south;
624  dsec = old * 10;
625  dsec2 = floor(dsec + 0.5);
626  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
627  s_off = diff;
628 
629  if (n_off < llepsilon || n_off <= s_off) {
630  old = cellhds.north;
631  dsec = old * 10;
632  dsec2 = floor(dsec + 0.5);
633  new = dsec2 / 10;
634  diff = n_off;
635  if (diff > 0 && diff < llepsilon) {
636  G_lat_format(old / 3600, buf);
637  G_lat_format(new / 3600, buf2);
638  if (strcmp(buf, buf2))
639  G_verbose_message(_("North rounded from %s to %s"),
640  buf, buf2);
641  cellhds.north = new;
642  }
643 
644  old = cellhds.south;
645  new = cellhds.north - cellhds.ns_res * cellhds.rows;
646  diff = fabs(new - old) / cellhds.ns_res;
647  if (diff > 0) {
648  G_lat_format(old / 3600, buf);
649  G_lat_format(new / 3600, buf2);
650  if (strcmp(buf, buf2))
651  G_verbose_message(_("South adjusted from %s to %s"),
652  buf, buf2);
653  }
654  cellhds.south = new;
655  }
656  else {
657  old = cellhds.south;
658  dsec = old * 10;
659  dsec2 = floor(dsec + 0.5);
660  new = dsec2 / 10;
661  diff = s_off;
662  if (diff > 0 && diff < llepsilon) {
663  G_lat_format(old / 3600, buf);
664  G_lat_format(new / 3600, buf2);
665  if (strcmp(buf, buf2))
666  G_verbose_message(_("South rounded from %s to %s"),
667  buf, buf2);
668  cellhds.south = new;
669  }
670 
671  old = cellhds.north;
672  new = cellhds.south + cellhds.ns_res * cellhds.rows;
673  diff = fabs(new - old) / cellhds.ns_res;
674  if (diff > 0) {
675  G_lat_format(old / 3600, buf);
676  G_lat_format(new / 3600, buf2);
677  if (strcmp(buf, buf2))
678  G_verbose_message(_("North adjusted from %s to %s"),
679  buf, buf2);
680  }
681  cellhds.north = new;
682  }
683  }
684  else {
685  old = cellhds.north;
686  dsec = old * 10;
687  dsec2 = floor(dsec + 0.5);
688  new = dsec2 / 10;
689  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
690  if (diff > 0 && diff < llepsilon) {
691  G_lat_format(old / 3600, buf);
692  G_lat_format(new / 3600, buf2);
693  if (strcmp(buf, buf2))
694  G_verbose_message(_("North rounded from %s to %s"),
695  buf, buf2);
696  ll_adjust = 1;
697  cellhds.north = new;
698  }
699 
700  old = cellhds.south;
701  dsec = old * 10;
702  dsec2 = floor(dsec + 0.5);
703  new = dsec2 / 10;
704  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
705  if (diff > 0 && diff < llepsilon) {
706  G_lat_format(old / 3600, buf);
707  G_lat_format(new / 3600, buf2);
708  if (strcmp(buf, buf2))
709  G_verbose_message(_("South rounded from %s to %s"),
710  buf, buf2);
711  ll_adjust = 1;
712  cellhds.south = new;
713  }
714  }
715  cellhds.ns_res = (cellhds.north - cellhds.south) / cellhds.rows;
716 
717  /* E - W */
718  /* resolution */
719  res_adj = 0;
720  old = cellhds.ew_res;
721 
722  if (old > 0.4) {
723  /* round to nearest 0.1 sec */
724  dsec = old * 10;
725  dsec2 = floor(dsec + 0.5);
726  new = dsec2 / 10;
727  diff = fabs(dsec2 - dsec) / dsec;
728  if (diff > 0 && diff < llepsilon) {
729  G_llres_format(old / 3600, buf);
730  G_llres_format(new / 3600, buf2);
731  if (strcmp(buf, buf2))
732  G_verbose_message(_("EW resolution rounded from %s to %s"),
733  buf, buf2);
734  ll_adjust = 1;
735  res_adj = 1;
736  cellhds.ew_res = new;
737  }
738  }
739 
740  if (res_adj) {
741  double w_off, e_off;
742 
743  old = cellhds.west;
744  dsec = old * 10;
745  dsec2 = floor(dsec + 0.5);
746  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
747  w_off = diff;
748 
749  old = cellhds.east;
750  dsec = old * 10;
751  dsec2 = floor(dsec + 0.5);
752  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
753  e_off = diff;
754 
755  if (w_off < llepsilon || w_off <= e_off) {
756  old = cellhds.west;
757  dsec = old * 10;
758  dsec2 = floor(dsec + 0.5);
759  new = dsec2 / 10;
760  diff = w_off;
761  if (diff > 0 && diff < llepsilon) {
762  G_lon_format(old / 3600, buf);
763  G_lon_format(new / 3600, buf2);
764  if (strcmp(buf, buf2))
765  G_verbose_message(_("West rounded from %s to %s"),
766  buf, buf2);
767  cellhds.west = new;
768  }
769 
770  old = cellhds.east;
771  new = cellhds.west + cellhds.ew_res * cellhds.cols;
772  diff = fabs(new - old) / cellhds.ew_res;
773  if (diff > 0) {
774  G_lon_format(old / 3600, buf);
775  G_lon_format(new / 3600, buf2);
776  if (strcmp(buf, buf2))
777  G_verbose_message(_("East adjusted from %s to %s"),
778  buf, buf2);
779  }
780  cellhds.east = new;
781  }
782  else {
783  old = cellhds.east;
784  dsec = old * 10;
785  dsec2 = floor(dsec + 0.5);
786  new = dsec2 / 10;
787  diff = e_off;
788  if (diff > 0 && diff < llepsilon) {
789  G_lon_format(old / 3600, buf);
790  G_lon_format(new / 3600, buf2);
791  if (strcmp(buf, buf2))
792  G_verbose_message(_("East rounded from %s to %s"),
793  buf, buf2);
794  cellhds.east = new;
795  }
796 
797  old = cellhds.west;
798  new = cellhds.east - cellhds.ew_res * cellhds.cols;
799  diff = fabs(new - cellhds.west) / cellhds.ew_res;
800  if (diff > 0) {
801  G_lon_format(old / 3600, buf);
802  G_lon_format(new / 3600, buf2);
803  if (strcmp(buf, buf2))
804  G_verbose_message(_("West adjusted from %s to %s"),
805  buf, buf2);
806  }
807  cellhds.west = new;
808  }
809  }
810  else {
811  old = cellhds.west;
812  dsec = old * 10;
813  dsec2 = floor(dsec + 0.5);
814  new = dsec2 / 10;
815  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
816  if (diff > 0 && diff < llepsilon) {
817  G_lon_format(old / 3600, buf);
818  G_lon_format(new / 3600, buf2);
819  if (strcmp(buf, buf2))
820  G_verbose_message(_("West rounded from %s to %s"),
821  buf, buf2);
822  ll_adjust = 1;
823  cellhds.west = new;
824  }
825 
826  old = cellhds.east;
827  dsec = old * 10;
828  dsec2 = floor(dsec + 0.5);
829  new = dsec2 / 10;
830  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
831  if (diff > 0 && diff < llepsilon) {
832  G_lon_format(old / 3600, buf);
833  G_lon_format(new / 3600, buf2);
834  if (strcmp(buf, buf2))
835  G_verbose_message(_("East rounded from %s to %s"),
836  buf, buf2);
837  ll_adjust = 1;
838  cellhds.east = new;
839  }
840  }
841  cellhds.ew_res = (cellhds.east - cellhds.west) / cellhds.cols;
842 
843  cellhd->ns_res = cellhds.ns_res / 3600;
844  cellhd->ew_res = cellhds.ew_res / 3600;
845  cellhd->north = cellhds.north / 3600;
846  cellhd->south = cellhds.south / 3600;
847  cellhd->west = cellhds.west / 3600;
848  cellhd->east = cellhds.east / 3600;
849 
850  return ll_adjust;
851 }
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
Definition: gis/error.c:130
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition: gis/error.c:108
if(!YY_CURRENT_BUFFER)
Definition: sqlp.yy.c:803
void G_adjust_Cell_head3(struct Cell_head *cellhd, int row_flag, int col_flag, int depth_flag)
Adjust cell header for 3D values.
Definition: adj_cellhd.c:150
2D/3D raster map header (used also for region)
Definition: gis.h:390
double west
Extent coordinates (west)
Definition: gis.h:442
void G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag)
Adjust cell header.
Definition: adj_cellhd.c:51
int old
Definition: raster3d/cats.c:84
int G_adjust_window_ll(struct Cell_head *cellhd)
Adjust window for lat/lon.
Definition: adj_cellhd.c:512
void G_lon_format(double lon, char *buf)
Definition: ll_format.c:56
#define LL_TOLERANCE
Definition: adj_cellhd.c:19
int cols3
Number of columns for 3D data.
Definition: gis.h:411
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
double top
Extent coordinates (top) - 3D data.
Definition: gis.h:444
double north
Extent coordinates (north)
Definition: gis.h:436
double ns_res3
Resolution - north to south cell size for 3D data.
Definition: gis.h:432
int rows3
Number of rows for 3D data.
Definition: gis.h:407
int G_lat_scan(const char *buf, double *lat)
Definition: ll_scan.c:46
double south
Extent coordinates (south)
Definition: gis.h:438
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double bottom
Extent coordinates (bottom) - 3D data.
Definition: gis.h:446
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:99
int depths
number of depths for 3D data
Definition: gis.h:413
int proj
Projection code.
Definition: gis.h:422
void G_llres_format(double res, char *buf)
Definition: ll_format.c:71
int cols
Number of columns for 2D data.
Definition: gis.h:409
int G_llres_scan(const char *buf, double *res)
Definition: ll_scan.c:56
double ns_res
Resolution - north to south cell size for 2D data.
Definition: gis.h:430
double east
Extent coordinates (east)
Definition: gis.h:440
#define _(str)
Definition: glocale.h:13
void G_lat_format(double lat, char *buf)
Definition: ll_format.c:41
double ew_res
Resolution - east to west cell size for 2D data.
Definition: gis.h:426
double tb_res
Resolution - top to bottom cell size for 3D data.
Definition: gis.h:434
int rows
Number of rows for 2D data.
Definition: gis.h:405
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
int G_lon_scan(const char *buf, double *lon)
Definition: ll_scan.c:51
double ew_res3
Resolution - east to west cell size for 3D data.
Definition: gis.h:428