16 #include <grass/gis.h>
17 #include <grass/glocale.h>
19 #define LL_TOLERANCE 10
23 static double llepsilon = 0.01;
24 static double fpepsilon = 1.0e-9;
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);
56 if (cellhd->ns_res <= 0)
57 G_fatal_error(_(
"Illegal n-s resolution value <%lf>"), cellhd->ns_res);
60 if (cellhd->rows <= 0)
64 if (cellhd->ew_res <= 0)
68 if (cellhd->cols <= 0)
73 if (cellhd->north <= cellhd->south) {
74 if (cellhd->proj == PROJECTION_LL)
82 if (cellhd->east <= cellhd->west)
88 (cellhd->north - cellhd->south +
89 cellhd->ns_res / 2.0) / cellhd->ns_res;
90 if (cellhd->rows == 0)
95 (cellhd->east - cellhd->west +
96 cellhd->ew_res / 2.0) / cellhd->ew_res;
97 if (cellhd->cols == 0)
101 if (cellhd->cols < 0 || cellhd->rows < 0) {
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)
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)
116 if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
151 int col_flag,
int depth_flag)
156 if (cellhd->ns_res <= 0)
158 if (cellhd->ns_res3 <= 0)
162 if (cellhd->rows <= 0)
164 if (cellhd->rows3 <= 0)
168 if (cellhd->ew_res <= 0)
170 if (cellhd->ew_res3 <= 0)
174 if (cellhd->cols <= 0)
176 if (cellhd->cols3 <= 0)
180 if (cellhd->tb_res <= 0)
184 if (cellhd->depths <= 0)
189 if (cellhd->north <= cellhd->south) {
190 if (cellhd->proj == PROJECTION_LL)
198 if (cellhd->east <= cellhd->west)
201 if (cellhd->top <= cellhd->bottom)
207 (cellhd->north - cellhd->south +
208 cellhd->ns_res / 2.0) / cellhd->ns_res;
209 if (cellhd->rows == 0)
213 (cellhd->north - cellhd->south +
214 cellhd->ns_res3 / 2.0) / cellhd->ns_res3;
215 if (cellhd->rows3 == 0)
220 (cellhd->east - cellhd->west +
221 cellhd->ew_res / 2.0) / cellhd->ew_res;
222 if (cellhd->cols == 0)
226 (cellhd->east - cellhd->west +
227 cellhd->ew_res3 / 2.0) / cellhd->ew_res3;
228 if (cellhd->cols3 == 0)
234 (cellhd->top - cellhd->bottom +
235 cellhd->tb_res / 2.0) / cellhd->tb_res;
236 if (cellhd->depths == 0)
240 if (cellhd->cols < 0 || cellhd->rows < 0 || cellhd->cols3 < 0 ||
241 cellhd->rows3 < 0 || cellhd->depths < 0) {
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)
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)
256 if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
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;
267 static int ll_wrap(
struct Cell_head *cellhd)
272 if (cellhd->proj != PROJECTION_LL)
275 if (cellhd->east <= cellhd->west) {
276 G_warning(_(
"East (%.15g) is not larger than West (%.15g)"),
277 cellhd->east, cellhd->west);
279 while (cellhd->east <= cellhd->west)
280 cellhd->east += 360.0;
289 while (cellhd->west + shift >= 180) {
292 while (cellhd->east + shift <= -180) {
297 while (cellhd->east + shift > 360) {
300 while (cellhd->west + shift <= -360) {
305 cellhd->west += shift;
306 cellhd->east += shift;
311 G_fatal_error(_(
"Illegal latitude for North: %g"), cellhd->north);
313 G_fatal_error(_(
"Illegal latitude for South: %g"), cellhd->south);
318 G_debug(1,
"East: %g", cellhd->east);
319 G_fatal_error(_(
"Illegal longitude for West: %g"), cellhd->west);
322 G_debug(1,
"West: %g", cellhd->west);
323 G_fatal_error(_(
"Illegal longitude for East: %g"), cellhd->east);
330 static int ll_check_ns(
struct Cell_head *cellhd)
337 if (cellhd->proj != PROJECTION_LL)
342 G_debug(3,
"ll_check_ns: epsilon: %g", llepsilon);
346 diff = (cellhd->north - cellhd->south) / cellhd->ns_res;
347 ncells = (
int) (diff + 0.5);
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"),
356 diff = (cellhd->north - 90) / cellhd->ns_res;
359 if (cellhd->north < 90.0 && diff < 1.0 ) {
362 if (diff < llepsilon && diff > fpepsilon) {
364 cellhd->north - 90.0);
371 if (cellhd->north > 90.0) {
372 if (diff <= 0.5 + llepsilon) {
376 if (diff < llepsilon && diff > fpepsilon) {
378 cellhd->north - 90.0);
379 G_debug(1,
"North of north in seconds: %g",
380 (cellhd->north - 90.0) * 3600);
390 if (diff < llepsilon && diff > fpepsilon) {
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);
406 diff = (cellhd->south + 90) / cellhd->ns_res;
409 if (cellhd->south > -90.0 && diff < 1.0 ) {
412 if (diff < llepsilon && diff > fpepsilon) {
414 cellhd->south + 90.0);
421 if (cellhd->south < -90.0) {
422 if (diff <= 0.5 + llepsilon) {
426 if (diff < llepsilon && diff > fpepsilon) {
429 G_debug(1,
"South of south in seconds: %g",
430 (-cellhd->south - 90) * 3600);
440 if (diff < llepsilon && diff > fpepsilon) {
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);
456 cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
461 static int ll_check_ew(
struct Cell_head *cellhd)
468 if (cellhd->proj != PROJECTION_LL)
473 G_debug(3,
"ll_check_ew: epsilon: %g", llepsilon);
476 diff = (cellhd->east - cellhd->west) / cellhd->ew_res;
477 ncells = (
int) (diff + 0.5);
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"),
484 if (cellhd->east - cellhd->west > 360.0) {
485 diff = (cellhd->east - cellhd->west - 360.0) / cellhd->ew_res;
486 if (diff > fpepsilon)
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)
514 int ll_adjust, res_adj;
516 char buf[100], buf2[100];
519 struct Cell_head cellhds;
522 if (cellhd->proj != PROJECTION_LL)
529 cellhd->ns_res =
new;
534 cellhd->ew_res =
new;
559 old = cellhds.ns_res * 3600;
560 sprintf(buf,
"%f", old);
561 sscanf(buf,
"%lf", &
new);
562 cellhds.ns_res =
new;
564 old = cellhds.ew_res * 3600;
565 sprintf(buf,
"%f", old);
566 sscanf(buf,
"%lf", &
new);
567 cellhds.ew_res =
new;
569 old = cellhds.north * 3600;
570 sprintf(buf,
"%f", old);
571 sscanf(buf,
"%lf", &
new);
574 old = cellhds.south * 3600;
575 sprintf(buf,
"%f", old);
576 sscanf(buf,
"%lf", &
new);
579 old = cellhds.west * 3600;
580 sprintf(buf,
"%f", old);
581 sscanf(buf,
"%lf", &
new);
584 old = cellhds.east * 3600;
585 sprintf(buf,
"%f", old);
586 sscanf(buf,
"%lf", &
new);
594 old = cellhds.ns_res;
599 dsec2 = floor(dsec + 0.5);
601 diff = fabs(dsec2 - dsec) / dsec;
602 if (diff > 0 && diff < llepsilon) {
605 if (strcmp(buf, buf2))
610 cellhds.ns_res =
new;
619 dsec2 = floor(dsec + 0.5);
620 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
625 dsec2 = floor(dsec + 0.5);
626 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
629 if (n_off < llepsilon || n_off <= s_off) {
632 dsec2 = floor(dsec + 0.5);
635 if (diff > 0 && diff < llepsilon) {
638 if (strcmp(buf, buf2))
645 new = cellhds.north - cellhds.ns_res * cellhds.rows;
646 diff = fabs(
new - old) / cellhds.ns_res;
650 if (strcmp(buf, buf2))
659 dsec2 = floor(dsec + 0.5);
662 if (diff > 0 && diff < llepsilon) {
665 if (strcmp(buf, buf2))
672 new = cellhds.south + cellhds.ns_res * cellhds.rows;
673 diff = fabs(
new - old) / cellhds.ns_res;
677 if (strcmp(buf, buf2))
687 dsec2 = floor(dsec + 0.5);
689 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
690 if (diff > 0 && diff < llepsilon) {
693 if (strcmp(buf, buf2))
702 dsec2 = floor(dsec + 0.5);
704 diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
705 if (diff > 0 && diff < llepsilon) {
708 if (strcmp(buf, buf2))
715 cellhds.ns_res = (cellhds.north - cellhds.south) / cellhds.rows;
720 old = cellhds.ew_res;
725 dsec2 = floor(dsec + 0.5);
727 diff = fabs(dsec2 - dsec) / dsec;
728 if (diff > 0 && diff < llepsilon) {
731 if (strcmp(buf, buf2))
736 cellhds.ew_res =
new;
745 dsec2 = floor(dsec + 0.5);
746 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
751 dsec2 = floor(dsec + 0.5);
752 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
755 if (w_off < llepsilon || w_off <= e_off) {
758 dsec2 = floor(dsec + 0.5);
761 if (diff > 0 && diff < llepsilon) {
764 if (strcmp(buf, buf2))
771 new = cellhds.west + cellhds.ew_res * cellhds.cols;
772 diff = fabs(
new - old) / cellhds.ew_res;
776 if (strcmp(buf, buf2))
785 dsec2 = floor(dsec + 0.5);
788 if (diff > 0 && diff < llepsilon) {
791 if (strcmp(buf, buf2))
798 new = cellhds.east - cellhds.ew_res * cellhds.cols;
799 diff = fabs(
new - cellhds.west) / cellhds.ew_res;
803 if (strcmp(buf, buf2))
813 dsec2 = floor(dsec + 0.5);
815 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
816 if (diff > 0 && diff < llepsilon) {
819 if (strcmp(buf, buf2))
828 dsec2 = floor(dsec + 0.5);
830 diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
831 if (diff > 0 && diff < llepsilon) {
834 if (strcmp(buf, buf2))
841 cellhds.ew_res = (cellhds.east - cellhds.west) / cellhds.cols;
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;