GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e384be2c9
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"));
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"));
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  diff);
530  }
531  else if (cellhd->east - cellhd->west < 360.0) {
532  diff = (360.0 - (cellhd->east - cellhd->west)) / cellhd->ew_res;
533  if (diff < 1.0 && diff > fpepsilon)
534  G_verbose_message(_("%g cells missing to cover 360 degree EW extent"),
535  diff);
536  }
537 
538  return lladjust;
539 }
540 
541 /*!
542  * \brief Adjust window for lat/lon.
543  *
544  * This function tries to automatically fix fp precision issues and
545  * adjust rounding errors for lat/lon.
546  *
547  * <b>Note:</b> 3D values are not adjusted.
548  *
549  * \param[in,out] cellhd pointer to Cell_head structure
550  * \return 1 if window was adjusted
551  * \return 0 if window was not adjusted
552  */
553 int G_adjust_window_ll(struct Cell_head *cellhd)
554 {
555  int ll_adjust, res_adj;
556  double dsec, dsec2;
557  char buf[100], buf2[100];
558  double diff;
559  double old, new;
560  struct Cell_head cellhds; /* everything in seconds, not degrees */
561 
562  /* lat/lon checks */
563  if (cellhd->proj != PROJECTION_LL)
564  return 0;
565 
566  /* put everything through ll_format + ll_scan */
567  G_llres_format(cellhd->ns_res, buf);
568  if (G_llres_scan(buf, &new) != 1)
569  G_fatal_error(_("Invalid NS resolution"));
570  cellhd->ns_res = new;
571 
572  G_llres_format(cellhd->ew_res, buf);
573  if (G_llres_scan(buf, &new) != 1)
574  G_fatal_error(_("Invalid EW resolution"));
575  cellhd->ew_res = new;
576 
577  G_lat_format(cellhd->north, buf);
578  if (G_lat_scan(buf, &new) != 1)
579  G_fatal_error(_("Invalid North"));
580  cellhd->north = new;
581 
582  G_lat_format(cellhd->south, buf);
583  if (G_lat_scan(buf, &new) != 1)
584  G_fatal_error(_("Invalid South"));
585  cellhd->south = new;
586 
587  G_lon_format(cellhd->west, buf);
588  if (G_lon_scan(buf, &new) != 1)
589  G_fatal_error(_("Invalid West"));
590  cellhd->west = new;
591 
592  G_lon_format(cellhd->east, buf);
593  if (G_lon_scan(buf, &new) != 1)
594  G_fatal_error(_("Invalid East"));
595  cellhd->east = new;
596 
597  /* convert to seconds */
598  cellhds = *cellhd;
599 
600  old = cellhds.ns_res * 3600;
601  sprintf(buf, "%f", old);
602  sscanf(buf, "%lf", &new);
603  cellhds.ns_res = new;
604 
605  old = cellhds.ew_res * 3600;
606  sprintf(buf, "%f", old);
607  sscanf(buf, "%lf", &new);
608  cellhds.ew_res = new;
609 
610  old = cellhds.north * 3600;
611  sprintf(buf, "%f", old);
612  sscanf(buf, "%lf", &new);
613  cellhds.north = new;
614 
615  old = cellhds.south * 3600;
616  sprintf(buf, "%f", old);
617  sscanf(buf, "%lf", &new);
618  cellhds.south = new;
619 
620  old = cellhds.west * 3600;
621  sprintf(buf, "%f", old);
622  sscanf(buf, "%lf", &new);
623  cellhds.west = new;
624 
625  old = cellhds.east * 3600;
626  sprintf(buf, "%f", old);
627  sscanf(buf, "%lf", &new);
628  cellhds.east = new;
629 
630  ll_adjust = 0;
631 
632  /* N - S */
633  /* resolution */
634  res_adj = 0;
635  old = cellhds.ns_res;
636 
637  if (old > 0.4) {
638  /* round to nearest 0.1 sec */
639  dsec = old * 10;
640  dsec2 = floor(dsec + 0.5);
641  new = dsec2 / 10;
642  diff = fabs(dsec2 - dsec) / dsec;
643  if (diff > 0 && diff < llepsilon) {
644  G_llres_format(old / 3600, buf);
645  G_llres_format(new / 3600, buf2);
646  if (strcmp(buf, buf2))
647  G_verbose_message(_("NS resolution rounded from %s to %s"),
648  buf, buf2);
649  ll_adjust = 1;
650  res_adj = 1;
651  cellhds.ns_res = new;
652  }
653  }
654 
655  if (res_adj) {
656  double n_off, s_off;
657 
658  old = cellhds.north;
659  dsec = old * 10;
660  dsec2 = floor(dsec + 0.5);
661  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
662  n_off = diff;
663 
664  old = cellhds.south;
665  dsec = old * 10;
666  dsec2 = floor(dsec + 0.5);
667  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
668  s_off = diff;
669 
670  if (n_off < llepsilon || n_off <= s_off) {
671  old = cellhds.north;
672  dsec = old * 10;
673  dsec2 = floor(dsec + 0.5);
674  new = dsec2 / 10;
675  diff = n_off;
676  if (diff > 0 && diff < llepsilon) {
677  G_lat_format(old / 3600, buf);
678  G_lat_format(new / 3600, buf2);
679  if (strcmp(buf, buf2))
680  G_verbose_message(_("North rounded from %s to %s"),
681  buf, buf2);
682  cellhds.north = new;
683  }
684 
685  old = cellhds.south;
686  new = cellhds.north - cellhds.ns_res * cellhds.rows;
687  diff = fabs(new - old) / cellhds.ns_res;
688  if (diff > 0) {
689  G_lat_format(old / 3600, buf);
690  G_lat_format(new / 3600, buf2);
691  if (strcmp(buf, buf2))
692  G_verbose_message(_("South adjusted from %s to %s"),
693  buf, buf2);
694  }
695  cellhds.south = new;
696  }
697  else {
698  old = cellhds.south;
699  dsec = old * 10;
700  dsec2 = floor(dsec + 0.5);
701  new = dsec2 / 10;
702  diff = s_off;
703  if (diff > 0 && diff < llepsilon) {
704  G_lat_format(old / 3600, buf);
705  G_lat_format(new / 3600, buf2);
706  if (strcmp(buf, buf2))
707  G_verbose_message(_("South rounded from %s to %s"),
708  buf, buf2);
709  cellhds.south = new;
710  }
711 
712  old = cellhds.north;
713  new = cellhds.south + cellhds.ns_res * cellhds.rows;
714  diff = fabs(new - old) / cellhds.ns_res;
715  if (diff > 0) {
716  G_lat_format(old / 3600, buf);
717  G_lat_format(new / 3600, buf2);
718  if (strcmp(buf, buf2))
719  G_verbose_message(_("North adjusted from %s to %s"),
720  buf, buf2);
721  }
722  cellhds.north = new;
723  }
724  }
725  else {
726  old = cellhds.north;
727  dsec = old * 10;
728  dsec2 = floor(dsec + 0.5);
729  new = dsec2 / 10;
730  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
731  if (diff > 0 && diff < llepsilon) {
732  G_lat_format(old / 3600, buf);
733  G_lat_format(new / 3600, buf2);
734  if (strcmp(buf, buf2))
735  G_verbose_message(_("North rounded from %s to %s"),
736  buf, buf2);
737  ll_adjust = 1;
738  cellhds.north = new;
739  }
740 
741  old = cellhds.south;
742  dsec = old * 10;
743  dsec2 = floor(dsec + 0.5);
744  new = dsec2 / 10;
745  diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
746  if (diff > 0 && diff < llepsilon) {
747  G_lat_format(old / 3600, buf);
748  G_lat_format(new / 3600, buf2);
749  if (strcmp(buf, buf2))
750  G_verbose_message(_("South rounded from %s to %s"),
751  buf, buf2);
752  ll_adjust = 1;
753  cellhds.south = new;
754  }
755  }
756  cellhds.ns_res = (cellhds.north - cellhds.south) / cellhds.rows;
757 
758  /* E - W */
759  /* resolution */
760  res_adj = 0;
761  old = cellhds.ew_res;
762 
763  if (old > 0.4) {
764  /* round to nearest 0.1 sec */
765  dsec = old * 10;
766  dsec2 = floor(dsec + 0.5);
767  new = dsec2 / 10;
768  diff = fabs(dsec2 - dsec) / dsec;
769  if (diff > 0 && diff < llepsilon) {
770  G_llres_format(old / 3600, buf);
771  G_llres_format(new / 3600, buf2);
772  if (strcmp(buf, buf2))
773  G_verbose_message(_("EW resolution rounded from %s to %s"),
774  buf, buf2);
775  ll_adjust = 1;
776  res_adj = 1;
777  cellhds.ew_res = new;
778  }
779  }
780 
781  if (res_adj) {
782  double w_off, e_off;
783 
784  old = cellhds.west;
785  dsec = old * 10;
786  dsec2 = floor(dsec + 0.5);
787  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
788  w_off = diff;
789 
790  old = cellhds.east;
791  dsec = old * 10;
792  dsec2 = floor(dsec + 0.5);
793  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
794  e_off = diff;
795 
796  if (w_off < llepsilon || w_off <= e_off) {
797  old = cellhds.west;
798  dsec = old * 10;
799  dsec2 = floor(dsec + 0.5);
800  new = dsec2 / 10;
801  diff = w_off;
802  if (diff > 0 && diff < llepsilon) {
803  G_lon_format(old / 3600, buf);
804  G_lon_format(new / 3600, buf2);
805  if (strcmp(buf, buf2))
806  G_verbose_message(_("West rounded from %s to %s"),
807  buf, buf2);
808  cellhds.west = new;
809  }
810 
811  old = cellhds.east;
812  new = cellhds.west + cellhds.ew_res * cellhds.cols;
813  diff = fabs(new - old) / cellhds.ew_res;
814  if (diff > 0) {
815  G_lon_format(old / 3600, buf);
816  G_lon_format(new / 3600, buf2);
817  if (strcmp(buf, buf2))
818  G_verbose_message(_("East adjusted from %s to %s"),
819  buf, buf2);
820  }
821  cellhds.east = new;
822  }
823  else {
824  old = cellhds.east;
825  dsec = old * 10;
826  dsec2 = floor(dsec + 0.5);
827  new = dsec2 / 10;
828  diff = e_off;
829  if (diff > 0 && diff < llepsilon) {
830  G_lon_format(old / 3600, buf);
831  G_lon_format(new / 3600, buf2);
832  if (strcmp(buf, buf2))
833  G_verbose_message(_("East rounded from %s to %s"),
834  buf, buf2);
835  cellhds.east = new;
836  }
837 
838  old = cellhds.west;
839  new = cellhds.east - cellhds.ew_res * cellhds.cols;
840  diff = fabs(new - cellhds.west) / cellhds.ew_res;
841  if (diff > 0) {
842  G_lon_format(old / 3600, buf);
843  G_lon_format(new / 3600, buf2);
844  if (strcmp(buf, buf2))
845  G_verbose_message(_("West adjusted from %s to %s"),
846  buf, buf2);
847  }
848  cellhds.west = new;
849  }
850  }
851  else {
852  old = cellhds.west;
853  dsec = old * 10;
854  dsec2 = floor(dsec + 0.5);
855  new = dsec2 / 10;
856  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
857  if (diff > 0 && diff < llepsilon) {
858  G_lon_format(old / 3600, buf);
859  G_lon_format(new / 3600, buf2);
860  if (strcmp(buf, buf2))
861  G_verbose_message(_("West rounded from %s to %s"),
862  buf, buf2);
863  ll_adjust = 1;
864  cellhds.west = new;
865  }
866 
867  old = cellhds.east;
868  dsec = old * 10;
869  dsec2 = floor(dsec + 0.5);
870  new = dsec2 / 10;
871  diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
872  if (diff > 0 && diff < llepsilon) {
873  G_lon_format(old / 3600, buf);
874  G_lon_format(new / 3600, buf2);
875  if (strcmp(buf, buf2))
876  G_verbose_message(_("East rounded from %s to %s"),
877  buf, buf2);
878  ll_adjust = 1;
879  cellhds.east = new;
880  }
881  }
882  cellhds.ew_res = (cellhds.east - cellhds.west) / cellhds.cols;
883 
884  cellhd->ns_res = cellhds.ns_res / 3600;
885  cellhd->ew_res = cellhds.ew_res / 3600;
886  cellhd->north = cellhds.north / 3600;
887  cellhd->south = cellhds.south / 3600;
888  cellhd->west = cellhds.west / 3600;
889  cellhd->east = cellhds.east / 3600;
890 
891  return ll_adjust;
892 }
if(!(yy_init))
Definition: sqlp.yy.c:775
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
Definition: gis/error.c:131
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition: gis/error.c:109
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:409
double west
Extent coordinates (west)
Definition: gis.h:461
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:553
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:430
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
double top
Extent coordinates (top) - 3D data.
Definition: gis.h:463
double north
Extent coordinates (north)
Definition: gis.h:455
double ns_res3
Resolution - north to south cell size for 3D data.
Definition: gis.h:451
int rows3
Number of rows for 3D data.
Definition: gis.h:426
int G_lat_scan(const char *buf, double *lat)
Definition: ll_scan.c:46
double south
Extent coordinates (south)
Definition: gis.h:457
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:465
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:113
int depths
number of depths for 3D data
Definition: gis.h:432
int proj
Projection code.
Definition: gis.h:441
void G_llres_format(double res, char *buf)
Definition: ll_format.c:71
int cols
Number of columns for 2D data.
Definition: gis.h:428
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:449
double east
Extent coordinates (east)
Definition: gis.h:459
#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:445
double tb_res
Resolution - top to bottom cell size for 3D data.
Definition: gis.h:453
int rows
Number of rows for 2D data.
Definition: gis.h:424
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
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:447