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