17#include <grass/glocale.h>
19#define LL_TOLERANCE 10
23static double llepsilon = 0.01;
24static double fpepsilon = 1.0e-9;
26static int ll_wrap(
struct Cell_head *cellhd);
27static int ll_check_ns(
struct Cell_head *cellhd);
28static int ll_check_ew(
struct Cell_head *cellhd);
56 if (cellhd->ns_res <= 0)
61 if (cellhd->rows <= 0)
63 " (resolution is %g)"),
64 cellhd->rows, cellhd->ns_res);
67 if (cellhd->ew_res <= 0)
72 if (cellhd->cols <= 0)
74 " (resolution is %g)"),
75 cellhd->cols, cellhd->ew_res);
79 if (cellhd->north <= cellhd->south) {
80 if (cellhd->proj == PROJECTION_LL)
82 " but %g (north) <= %g (south"),
83 cellhd->north, cellhd->south);
86 " but %g (north) <= %g (south"),
87 cellhd->north, cellhd->south);
92 if (cellhd->east <= cellhd->west)
94 " but %g (east) <= %g (west)"),
95 cellhd->east, cellhd->west);
99 cellhd->rows = (cellhd->north - cellhd->south + cellhd->ns_res / 2.0) /
101 if (cellhd->rows == 0)
105 cellhd->cols = (cellhd->east - cellhd->west + cellhd->ew_res / 2.0) /
107 if (cellhd->cols == 0)
111 if (cellhd->cols < 0) {
112 G_fatal_error(_(
"Invalid coordinates: negative number of columns"));
114 if (cellhd->rows < 0) {
115 G_fatal_error(_(
"Invalid coordinates: negative number of rows"));
119 old_res = cellhd->ns_res;
120 cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
121 if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
124 old_res = cellhd->ew_res;
125 cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
126 if (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
129 if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
169 if (cellhd->ns_res <= 0)
172 if (cellhd->ns_res3 <= 0)
177 if (cellhd->rows <= 0)
179 " (resolution is %g)"),
180 cellhd->rows, cellhd->ns_res);
181 if (cellhd->rows3 <= 0)
183 " (resolution is %g)"),
184 cellhd->rows3, cellhd->ns_res3);
187 if (cellhd->ew_res <= 0)
190 if (cellhd->ew_res3 <= 0)
195 if (cellhd->cols <= 0)
197 " (resolution is %g)"),
198 cellhd->cols, cellhd->ew_res);
199 if (cellhd->cols3 <= 0)
201 " (resolution is %g)"),
202 cellhd->cols3, cellhd->ew_res3);
205 if (cellhd->tb_res <= 0)
210 if (cellhd->depths <= 0)
211 G_fatal_error(_(
"Illegal depths value: %d"), cellhd->depths);
215 if (cellhd->north <= cellhd->south) {
216 if (cellhd->proj == PROJECTION_LL)
218 " but %g (north) <= %g (south"),
219 cellhd->north, cellhd->south);
222 " but %g (north) <= %g (south"),
223 cellhd->north, cellhd->south);
228 if (cellhd->east <= cellhd->west)
230 " but %g (east) <= %g (west)"),
231 cellhd->east, cellhd->west);
233 if (cellhd->top <= cellhd->bottom)
235 " but %g (top) <= %g (bottom)"),
236 cellhd->top, cellhd->bottom);
240 cellhd->rows = (cellhd->north - cellhd->south + cellhd->ns_res / 2.0) /
242 if (cellhd->rows == 0)
246 (cellhd->north - cellhd->south + cellhd->ns_res3 / 2.0) /
248 if (cellhd->rows3 == 0)
252 cellhd->cols = (cellhd->east - cellhd->west + cellhd->ew_res / 2.0) /
254 if (cellhd->cols == 0)
257 cellhd->cols3 = (cellhd->east - cellhd->west + cellhd->ew_res3 / 2.0) /
259 if (cellhd->cols3 == 0)
264 cellhd->depths = (cellhd->top - cellhd->bottom + cellhd->tb_res / 2.0) /
266 if (cellhd->depths == 0)
270 if (cellhd->cols < 0 || cellhd->cols3 < 0) {
271 G_fatal_error(_(
"Invalid coordinates: negative number of columns"));
273 if (cellhd->rows < 0 || cellhd->rows3 < 0) {
274 G_fatal_error(_(
"Invalid coordinates: negative number of rows"));
276 if (cellhd->depths < 0) {
277 G_fatal_error(_(
"Invalid coordinates: negative number of depths"));
281 old_res = cellhd->ns_res;
282 cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
283 if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
286 old_res = cellhd->ew_res;
287 cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
288 if (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
291 if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
297 cellhd->ns_res3 = (cellhd->north - cellhd->south) / cellhd->rows3;
298 cellhd->ew_res3 = (cellhd->east - cellhd->west) / cellhd->cols3;
299 cellhd->tb_res = (cellhd->top - cellhd->bottom) / cellhd->depths;
302static int ll_wrap(
struct Cell_head *cellhd)
307 if (cellhd->proj != PROJECTION_LL)
310 if (cellhd->east <= cellhd->west) {
311 G_warning(_(
"East (%.15g) is not larger than West (%.15g)"),
312 cellhd->east, cellhd->west);
314 while (cellhd->east <= cellhd->west)
315 cellhd->east += 360.0;
324 while (cellhd->west + shift >= 180) {
327 while (cellhd->east + shift <= -180) {
332 while (cellhd->east + shift > 360) {
335 while (cellhd->west + shift <= -360) {
340 cellhd->west += shift;
341 cellhd->east += shift;
346 G_fatal_error(_(
"Illegal latitude for North: %g"), cellhd->north);
348 G_fatal_error(_(
"Illegal latitude for South: %g"), cellhd->south);
353 G_debug(1,
"East: %g", cellhd->east);
354 G_fatal_error(_(
"Illegal longitude for West: %g"), cellhd->west);
357 G_debug(1,
"West: %g", cellhd->west);
358 G_fatal_error(_(
"Illegal longitude for East: %g"), cellhd->east);
365static int ll_check_ns(
struct Cell_head *cellhd)
372 if (cellhd->proj != PROJECTION_LL)
377 G_debug(3,
"ll_check_ns: epsilon: %g", llepsilon);
381 diff = (cellhd->north - cellhd->south) / cellhd->ns_res;
382 ncells = (int)(diff + 0.5);
384 if ((diff < 0 && diff < -fpepsilon) || (diff > 0 && diff > fpepsilon)) {
386 _(
"NS extent does not match NS resolution: %g cells difference"),
391 diff = (cellhd->north - 90) / cellhd->ns_res;
394 if (cellhd->north < 90.0 && diff < 1.0) {
396 if (diff < llepsilon && diff > fpepsilon) {
398 _(
"Subtle input data rounding error of north boundary (%g)"),
399 cellhd->north - 90.0);
406 if (cellhd->north > 90.0) {
407 if (diff <= 0.5 + llepsilon) {
411 if (diff < llepsilon && diff > fpepsilon) {
414 cellhd->north - 90.0);
415 G_debug(1,
"North of north in seconds: %g",
416 (cellhd->north - 90.0) * 3600);
426 if (diff < llepsilon && diff > fpepsilon) {
429 cellhd->north - 90.0 - cellhd->ns_res / 2.0);
430 G_debug(1,
"North of north + 0.5 cells in seconds: %g",
431 (cellhd->north - 90.0 - cellhd->ns_res / 2.0) * 3600);
439 G_fatal_error(_(
"Illegal latitude for North: %g"), cellhd->north);
443 diff = (cellhd->south + 90) / cellhd->ns_res;
446 if (cellhd->south > -90.0 && diff < 1.0) {
448 if (diff < llepsilon && diff > fpepsilon) {
450 _(
"Subtle input data rounding error of south boundary (%g)"),
451 cellhd->south + 90.0);
458 if (cellhd->south < -90.0) {
459 if (diff <= 0.5 + llepsilon) {
463 if (diff < llepsilon && diff > fpepsilon) {
467 G_debug(1,
"South of south in seconds: %g",
468 (-cellhd->south - 90) * 3600);
478 if (diff < llepsilon && diff > fpepsilon) {
481 cellhd->south + 90 + cellhd->ns_res / 2.0);
482 G_debug(1,
"South of south + 0.5 cells in seconds: %g",
483 (-cellhd->south - 90 - cellhd->ns_res / 2.0) * 3600);
491 G_fatal_error(_(
"Illegal latitude for South: %g"), cellhd->south);
495 cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
500static int ll_check_ew(
struct Cell_head *cellhd)
507 if (cellhd->proj != PROJECTION_LL)
512 G_debug(3,
"ll_check_ew: epsilon: %g", llepsilon);
515 diff = (cellhd->east - cellhd->west) / cellhd->ew_res;
516 ncells = (int)(diff + 0.5);
518 if ((diff < 0 && diff < -fpepsilon) || (diff > 0 && diff > fpepsilon)) {
520 _(
"EW extent does not match EW resolution: %g cells difference"),
523 if (cellhd->east - cellhd->west > 360.0) {
524 diff = (cellhd->east - cellhd->west - 360.0) / cellhd->ew_res;
525 if (diff > fpepsilon)
527 " (East: %g, West: %g)"),
528 diff, cellhd->east, cellhd->west);
530 else if (cellhd->east - cellhd->west < 360.0) {
531 diff = (360.0 - (cellhd->east - cellhd->west)) / cellhd->ew_res;
532 if (diff < 1.0 && diff > fpepsilon)
534 _(
"%g cells missing to cover 360 degree EW extent"), diff);
554 int ll_adjust, res_adj;
556 char buf[100], buf2[100];
559 struct Cell_head cellhds;
562 if (cellhd->proj != PROJECTION_LL)
569 cellhd->ns_res =
new;
574 cellhd->ew_res =
new;
599 old = cellhds.ns_res * 3600;
600 sprintf(buf,
"%f", old);
601 sscanf(buf,
"%lf", &
new);
602 cellhds.ns_res =
new;
604 old = cellhds.ew_res * 3600;
605 sprintf(buf,
"%f", old);
606 sscanf(buf,
"%lf", &
new);
607 cellhds.ew_res =
new;
609 old = cellhds.north * 3600;
610 sprintf(buf,
"%f", old);
611 sscanf(buf,
"%lf", &
new);
614 old = cellhds.south * 3600;
615 sprintf(buf,
"%f", old);
616 sscanf(buf,
"%lf", &
new);
619 old = cellhds.west * 3600;
620 sprintf(buf,
"%f", old);
621 sscanf(buf,
"%lf", &
new);
624 old = cellhds.east * 3600;
625 sprintf(buf,
"%f", old);
626 sscanf(buf,
"%lf", &
new);
634 old = cellhds.ns_res;
639 dsec2 = floor(dsec + 0.5);
641 diff = fabs(dsec2 - dsec) / dsec;
642 if (diff > 0 && diff < llepsilon) {
645 if (strcmp(buf, buf2))
650 cellhds.ns_res =
new;
659 dsec2 = floor(dsec + 0.5);
660 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
665 dsec2 = floor(dsec + 0.5);
666 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
669 if (n_off < llepsilon || n_off <= s_off) {
672 dsec2 = floor(dsec + 0.5);
675 if (diff > 0 && diff < llepsilon) {
678 if (strcmp(buf, buf2))
685 new = cellhds.north - cellhds.ns_res *cellhds.rows;
686 diff = fabs(
new - old) / cellhds.ns_res;
690 if (strcmp(buf, buf2))
699 dsec2 = floor(dsec + 0.5);
702 if (diff > 0 && diff < llepsilon) {
705 if (strcmp(buf, buf2))
712 new = cellhds.south + cellhds.ns_res *cellhds.rows;
713 diff = fabs(
new - old) / cellhds.ns_res;
717 if (strcmp(buf, buf2))
727 dsec2 = floor(dsec + 0.5);
729 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
730 if (diff > 0 && diff < llepsilon) {
733 if (strcmp(buf, buf2))
741 dsec2 = floor(dsec + 0.5);
743 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
744 if (diff > 0 && diff < llepsilon) {
747 if (strcmp(buf, buf2))
753 cellhds.ns_res = (cellhds.north - cellhds.south) / cellhds.rows;
758 old = cellhds.ew_res;
763 dsec2 = floor(dsec + 0.5);
765 diff = fabs(dsec2 - dsec) / dsec;
766 if (diff > 0 && diff < llepsilon) {
769 if (strcmp(buf, buf2))
774 cellhds.ew_res =
new;
783 dsec2 = floor(dsec + 0.5);
784 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
789 dsec2 = floor(dsec + 0.5);
790 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
793 if (w_off < llepsilon || w_off <= e_off) {
796 dsec2 = floor(dsec + 0.5);
799 if (diff > 0 && diff < llepsilon) {
802 if (strcmp(buf, buf2))
809 new = cellhds.west + cellhds.ew_res *cellhds.cols;
810 diff = fabs(
new - old) / cellhds.ew_res;
814 if (strcmp(buf, buf2))
823 dsec2 = floor(dsec + 0.5);
826 if (diff > 0 && diff < llepsilon) {
829 if (strcmp(buf, buf2))
836 new = cellhds.east - cellhds.ew_res *cellhds.cols;
837 diff = fabs(
new - cellhds.west) / cellhds.ew_res;
841 if (strcmp(buf, buf2))
851 dsec2 = floor(dsec + 0.5);
853 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
854 if (diff > 0 && diff < llepsilon) {
857 if (strcmp(buf, buf2))
865 dsec2 = floor(dsec + 0.5);
867 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
868 if (diff > 0 && diff < llepsilon) {
871 if (strcmp(buf, buf2))
877 cellhds.ew_res = (cellhds.east - cellhds.west) / cellhds.cols;
879 cellhd->ns_res = cellhds.ns_res / 3600;
880 cellhd->ew_res = cellhds.ew_res / 3600;
881 cellhd->north = cellhds.north / 3600;
882 cellhd->south = cellhds.south / 3600;
883 cellhd->west = cellhds.west / 3600;
884 cellhd->east = cellhds.east / 3600;
void G_adjust_Cell_head3(struct Cell_head *cellhd, int row_flag, int col_flag, int depth_flag)
Adjust cell header for 3D values.
void G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag)
Adjust cell header.
int G_adjust_window_ll(struct Cell_head *cellhd)
Adjust window for lat/lon.
int G_debug(int level, const char *msg,...)
Print debugging message.
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
void G_warning(const char *msg,...)
Print a warning message to stderr.
int G_lat_scan(const char *buf, double *lat)
int G_lon_scan(const char *buf, double *lon)
int G_llres_scan(const char *buf, double *res)