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 */
G_adjust_Cell_head(struct Cell_head * cellhd,int row_flag,int col_flag)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 */
G_adjust_Cell_head3(struct Cell_head * cellhd,int row_flag,int col_flag,int depth_flag)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
ll_wrap(struct Cell_head * cellhd)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
ll_check_ns(struct Cell_head * cellhd)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
ll_check_ew(struct Cell_head * cellhd)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 */
G_adjust_window_ll(struct Cell_head * cellhd)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; /* everything 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 }
852