[GRASS-SVN] r70645 - grass/trunk/lib/gis
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Feb 20 06:18:21 PST 2017
Author: mmetz
Date: 2017-02-20 06:18:21 -0800 (Mon, 20 Feb 2017)
New Revision: 70645
Modified:
grass/trunk/lib/gis/adj_cellhd.c
Log:
libgis: improve G_adjust_Cell_head[3]() for latlong
Modified: grass/trunk/lib/gis/adj_cellhd.c
===================================================================
--- grass/trunk/lib/gis/adj_cellhd.c 2017-02-20 14:16:20 UTC (rev 70644)
+++ grass/trunk/lib/gis/adj_cellhd.c 2017-02-20 14:18:21 UTC (rev 70645)
@@ -11,9 +11,20 @@
* \author Original author CERL
*/
+#include <math.h>
+#include <string.h>
#include <grass/gis.h>
#include <grass/glocale.h>
+/* TODO: find good thresholds */
+/* deviation measured in cells */
+static double llepsilon = 0.01;
+static double fpepsilon = 1.0e-9;
+
+static int ll_wrap(struct Cell_head *cellhd);
+static int ll_check_ns(struct Cell_head *cellhd);
+static int ll_check_ew(struct Cell_head *cellhd);
+
/*!
* \brief Adjust cell header.
*
@@ -37,6 +48,8 @@
*/
void G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag)
{
+ double old_res;
+
if (!row_flag) {
if (cellhd->ns_res <= 0)
G_fatal_error(_("Illegal n-s resolution value <%lf>"), cellhd->ns_res);
@@ -54,72 +67,6 @@
G_fatal_error(_("Illegal col value"));
}
- /* for lat/lon, check north,south. force east larger than west */
- if (cellhd->proj == PROJECTION_LL) {
- double epsilon_ns, epsilon_ew;
-
- /* TODO: find good thresholds */
- epsilon_ns = 1. / cellhd->rows * 0.001;
- epsilon_ew = .000001; /* epsilon_ew calculation doesn't work due to cellhd->cols update/global wraparound below */
-
- G_debug(3, "G_adjust_Cell_head: epsilon_ns: %g, epsilon_ew: %g",
- epsilon_ns, epsilon_ew);
-
- /* TODO: once working, change below G_warning to G_debug */
-
- /* fix rounding problems if input map slightly exceeds the world definition -180 90 180 -90 */
- if (cellhd->north > 90.0) {
- if (((cellhd->north - 90.0) < epsilon_ns) &&
- ((cellhd->north - 90.0) > GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of north boundary (%g>%g)"),
- cellhd->north - 90.0, epsilon_ns);
- cellhd->north = 90.0;
- }
- else
- G_fatal_error(_("Illegal latitude for North"));
- }
-
- if (cellhd->south < -90.0) {
- if (((cellhd->south + 90.0) < epsilon_ns) &&
- ((cellhd->south + 90.0) < GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of south boundary (%g>%g)"),
- cellhd->south + 90.0, epsilon_ns);
- cellhd->south = -90.0;
- }
- else
- G_fatal_error(_("Illegal latitude for South"));
- }
-
-#if 0
- /* DISABLED: this breaks global wrap-around */
-
- G_debug(3,
- "G_adjust_Cell_head() cellhd->west: %f, devi: %g, eps: %g",
- cellhd->west, cellhd->west + 180.0, epsilon_ew);
-
- if ((cellhd->west < -180.0) && ((cellhd->west + 180.0) < epsilon_ew)
- && ((cellhd->west + 180.0) < GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of west boundary (%g>%g)"),
- cellhd->west + 180.0, epsilon_ew);
- cellhd->west = -180.0;
- }
-
- G_debug(3,
- "G_adjust_Cell_head() cellhd->east: %f, devi: %g, eps: %g",
- cellhd->east, cellhd->east - 180.0, epsilon_ew);
-
- if ((cellhd->east > 180.0) && ((cellhd->east - 180.0) > epsilon_ew)
- && ((cellhd->east - 180.0) > GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of east boundary (%g>%g)"),
- cellhd->east - 180.0, epsilon_ew);
- cellhd->east = 180.0;
- }
-#endif
-
- while (cellhd->east <= cellhd->west)
- cellhd->east += 360.0;
- }
-
/* check the edge values */
if (cellhd->north <= cellhd->south) {
if (cellhd->proj == PROJECTION_LL)
@@ -127,6 +74,9 @@
else
G_fatal_error(_("North must be larger than South"));
}
+
+ ll_wrap(cellhd);
+
if (cellhd->east <= cellhd->west)
G_fatal_error(_("East must be larger than West"));
@@ -150,10 +100,22 @@
G_fatal_error(_("Invalid coordinates"));
}
-
/* (re)compute the resolutions */
+ old_res = cellhd->ns_res;
cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
+ if (old_res > 0 && (old_res - cellhd->ns_res) / old_res > 0.001)
+ G_message(_("NS resolution has been changed"));
+
+ old_res = cellhd->ew_res;
cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
+ if (old_res > 0 && (old_res - cellhd->ew_res) / old_res > 0.001)
+ G_message(_("EW resolution has been changed"));
+
+ if ((cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.001)
+ G_message(_("NS and EW resolutions are different"));
+
+ ll_check_ns(cellhd);
+ ll_check_ew(cellhd);
}
/*!
@@ -186,6 +148,8 @@
void G_adjust_Cell_head3(struct Cell_head *cellhd, int row_flag,
int col_flag, int depth_flag)
{
+ double old_res;
+
if (!row_flag) {
if (cellhd->ns_res <= 0)
G_fatal_error(_("Illegal n-s resolution value"));
@@ -219,72 +183,6 @@
G_fatal_error(_("Illegal depths value"));
}
- /* for lat/lon, check north,south. force east larger than west */
- if (cellhd->proj == PROJECTION_LL) {
- double epsilon_ns, epsilon_ew;
-
- /* TODO: find good thresholds */
- epsilon_ns = 1. / cellhd->rows * 0.001;
- epsilon_ew = .000001; /* epsilon_ew calculation doesn't work due to cellhd->cols update/global wraparound below */
-
- G_debug(3, "G_adjust_Cell_head: epsilon_ns: %g, epsilon_ew: %g",
- epsilon_ns, epsilon_ew);
-
- /* TODO: once working, change below G_warning to G_debug */
-
- /* fix rounding problems if input map slightly exceeds the world definition -180 90 180 -90 */
- if (cellhd->north > 90.0) {
- if (((cellhd->north - 90.0) < epsilon_ns) &&
- ((cellhd->north - 90.0) > GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of north boundary (%g>%g)"),
- cellhd->north - 90.0, epsilon_ns);
- cellhd->north = 90.0;
- }
- else
- G_fatal_error(_("Illegal latitude for North"));
- }
-
- if (cellhd->south < -90.0) {
- if (((cellhd->south + 90.0) < epsilon_ns) &&
- ((cellhd->south + 90.0) < GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of south boundary (%g>%g)"),
- cellhd->south + 90.0, epsilon_ns);
- cellhd->south = -90.0;
- }
- else
- G_fatal_error(_("Illegal latitude for South"));
- }
-
-#if 0
- /* DISABLED: this breaks global wrap-around */
-
- G_debug(3,
- "G_adjust_Cell_head3() cellhd->west: %f, devi: %g, eps: %g",
- cellhd->west, cellhd->west + 180.0, epsilon_ew);
-
- if ((cellhd->west < -180.0) && ((cellhd->west + 180.0) < epsilon_ew)
- && ((cellhd->west + 180.0) < GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of west boundary (%g>%g)"),
- cellhd->west + 180.0, epsilon_ew);
- cellhd->west = -180.0;
- }
-
- G_debug(3,
- "G_adjust_Cell_head3() cellhd->east: %f, devi: %g, eps: %g",
- cellhd->east, cellhd->east - 180.0, epsilon_ew);
-
- if ((cellhd->east > 180.0) && ((cellhd->east - 180.0) > epsilon_ew)
- && ((cellhd->east - 180.0) > GRASS_EPSILON)) {
- G_warning(_("Fixing subtle input data rounding error of east boundary (%g>%g)"),
- cellhd->east - 180.0, epsilon_ew);
- cellhd->east = 180.0;
- }
-#endif
-
- while (cellhd->east <= cellhd->west)
- cellhd->east += 360.0;
- }
-
/* check the edge values */
if (cellhd->north <= cellhd->south) {
if (cellhd->proj == PROJECTION_LL)
@@ -292,6 +190,9 @@
else
G_fatal_error(_("North must be larger than South"));
}
+
+ ll_wrap(cellhd);
+
if (cellhd->east <= cellhd->west)
G_fatal_error(_("East must be larger than West"));
@@ -332,7 +233,6 @@
cellhd->tb_res / 2.0) / cellhd->tb_res;
if (cellhd->depths == 0)
cellhd->depths = 1;
-
}
if (cellhd->cols < 0 || cellhd->rows < 0 || cellhd->cols3 < 0 ||
@@ -341,9 +241,256 @@
}
/* (re)compute the resolutions */
+ old_res = cellhd->ns_res;
cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
+ if (old_res > 0 && (old_res - cellhd->ns_res) / old_res > 0.001)
+ G_message(_("NS resolution has been changed"));
+
+ old_res = cellhd->ew_res;
+ cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
+ if (old_res > 0 && (old_res - cellhd->ew_res) / old_res > 0.001)
+ G_message(_("EW resolution has been changed"));
+
+ if ((cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.001)
+ G_message(_("NS and EW resolutions are different"));
+
+ ll_check_ns(cellhd);
+ ll_check_ew(cellhd);
+
cellhd->ns_res3 = (cellhd->north - cellhd->south) / cellhd->rows3;
- cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
cellhd->ew_res3 = (cellhd->east - cellhd->west) / cellhd->cols3;
cellhd->tb_res = (cellhd->top - cellhd->bottom) / cellhd->depths;
}
+
+static int ll_wrap(struct Cell_head *cellhd)
+{
+ double shift;
+
+ /* for lat/lon, force east larger than west, try to wrap to -180, 180 */
+ if (cellhd->proj != PROJECTION_LL)
+ return 0;
+
+ if (cellhd->east <= cellhd->west) {
+ G_warning(_("East (%.15g) is not larger than West (%.15g)"),
+ cellhd->east, cellhd->west);
+
+ while (cellhd->east <= cellhd->west)
+ cellhd->east += 360.0;
+ }
+
+ /* with east larger than west,
+ * any 360 degree W-E extent can be represented within -360, 360
+ * but not within -180, 180 */
+
+ /* try to shift to within -180, 180 */
+ shift = 0;
+ while (cellhd->west + shift >= 180) {
+ shift -= 360.0;
+ }
+ while (cellhd->east + shift <= -180) {
+ shift += 360.0;
+ }
+
+ /* try to shift to within -360, 360 */
+ while (cellhd->east + shift > 360) {
+ shift -= 360.0;
+ }
+ while (cellhd->west + shift <= -360) {
+ shift += 360.0;
+ }
+
+ if (shift) {
+ cellhd->west += shift;
+ cellhd->east += shift;
+ }
+
+ /* very liberal thresholds */
+ if (cellhd->north > 100.0)
+ G_fatal_error(_("Illegal latitude for North: %g"), cellhd->north);
+ if (cellhd->south < -100.0)
+ G_fatal_error(_("Illegal latitude for South: %g"), cellhd->south);
+
+#if 0
+ /* disabled: allow W-E extents larger than 360 degree e.g. for display */
+ if (cellhd->west < -370.0) {
+ G_debug(1, "East: %g", cellhd->east);
+ G_fatal_error(_("Illegal longitude for West: %g"), cellhd->west);
+ }
+ if (cellhd->east > 370.0) {
+ G_debug(1, "West: %g", cellhd->west);
+ G_fatal_error(_("Illegal longitude for East: %g"), cellhd->east);
+ }
+#endif
+
+ return 1;
+}
+
+static int ll_check_ns(struct Cell_head *cellhd)
+{
+ int lladjust;
+ double diff;
+ int ncells;
+
+ /* lat/lon checks */
+ if (cellhd->proj != PROJECTION_LL)
+ return 0;
+
+ lladjust = 0;
+
+ G_debug(3, "ll_check_ns: epsilon: %g", llepsilon);
+
+ /* North, South: allow a half cell spill-over */
+
+ diff = (cellhd->north - cellhd->south) / cellhd->ns_res;
+ ncells = (int) (diff + 0.5);
+ diff -= ncells;
+ if ((diff < 0 && diff < -fpepsilon) ||
+ (diff > 0 && diff > fpepsilon)) {
+ G_message(_("NS extent does not match NS resolution: %g cells difference"),
+ diff);
+ }
+
+ /* north */
+ diff = (cellhd->north - 90) / cellhd->ns_res;
+ if (diff < 0)
+ diff = -diff;
+ if (cellhd->north < 90.0 && diff < 1.0 ) {
+ G_message(_("%g cells missing to reach 90 degree north"),
+ diff);
+ if (diff < llepsilon && diff > fpepsilon) {
+ G_message(_("Subtle input data rounding error of north boundary (%g)"),
+ cellhd->north - 90.0);
+ /* check only, do not modify
+ cellhd->north = 90.0;
+ lladjust = 1;
+ */
+ }
+ }
+ if (cellhd->north > 90.0) {
+ if (diff <= 0.5 + llepsilon) {
+ G_message(_("90 degree north is exceeded by %g cells"),
+ diff);
+
+ if (diff < llepsilon && diff > fpepsilon) {
+ G_message(_("Subtle input data rounding error of north boundary (%g)"),
+ cellhd->north - 90.0);
+ G_debug(1, "North of north in seconds: %g",
+ (cellhd->north - 90.0) * 3600);
+ /* check only, do not modify
+ cellhd->north = 90.0;
+ lladjust = 1;
+ */
+ }
+
+ diff = diff - 0.5;
+ if (diff < 0)
+ diff = -diff;
+ if (diff < llepsilon && diff > fpepsilon) {
+ G_message(_("Subtle input data rounding error of north boundary (%g)"),
+ cellhd->north - 90.0 - cellhd->ns_res / 2.0);
+ G_debug(1, "North of north + 0.5 cells in seconds: %g",
+ (cellhd->north - 90.0 - cellhd->ns_res / 2.0) * 3600);
+ /* check only, do not modify
+ cellhd->north = 90.0 + cellhd->ns_res / 2.0;
+ lladjust = 1;
+ */
+ }
+ }
+ else
+ G_fatal_error(_("Illegal latitude for North"));
+ }
+
+ /* south */
+ diff = (cellhd->south + 90) / cellhd->ns_res;
+ if (diff < 0)
+ diff = -diff;
+ if (cellhd->south > -90.0 && diff < 1.0 ) {
+ G_message(_("%g cells missing to reach 90 degree south"),
+ diff);
+ if (diff < llepsilon && diff > fpepsilon) {
+ G_message(_("Subtle input data rounding error of south boundary (%g)"),
+ cellhd->south + 90.0);
+ /* check only, do not modify
+ cellhd->south = -90.0;
+ lladjust = 1;
+ */
+ }
+ }
+ if (cellhd->south < -90.0) {
+ if (diff <= 0.5 + llepsilon) {
+ G_message(_("90 degree south is exceeded by %g cells"),
+ diff);
+
+ if (diff < llepsilon && diff > fpepsilon) {
+ G_message(_("Subtle input data rounding error of south boundary (%g)"),
+ cellhd->south + 90);
+ G_debug(1, "South of south in seconds: %g",
+ (-cellhd->south - 90) * 3600);
+ /* check only, do not modify
+ cellhd->south = -90.0;
+ lladjust = 1;
+ */
+ }
+
+ diff = diff - 0.5;
+ if (diff < 0)
+ diff = -diff;
+ if (diff < llepsilon && diff > fpepsilon) {
+ G_message(_("Subtle input data rounding error of south boundary (%g)"),
+ cellhd->south + 90 + cellhd->ns_res / 2.0);
+ G_debug(1, "South of south + 0.5 cells in seconds: %g",
+ (-cellhd->south - 90 - cellhd->ns_res / 2.0) * 3600);
+ /* check only, do not modify
+ cellhd->south = -90.0 - cellhd->ns_res / 2.0;
+ lladjust = 1;
+ */
+ }
+ }
+ else
+ G_fatal_error(_("Illegal latitude for South"));
+ }
+
+ if (lladjust)
+ cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
+
+ return lladjust;
+}
+
+static int ll_check_ew(struct Cell_head *cellhd)
+{
+ int lladjust;
+ double diff;
+ int ncells;
+
+ /* lat/lon checks */
+ if (cellhd->proj != PROJECTION_LL)
+ return 0;
+
+ lladjust = 0;
+
+ G_debug(3, "ll_check_ew: epsilon: %g", llepsilon);
+
+ /* west - east, no adjustment */
+ diff = (cellhd->east - cellhd->west) / cellhd->ew_res;
+ ncells = (int) (diff + 0.5);
+ diff -= ncells;
+ if ((diff < 0 && diff < -fpepsilon) ||
+ (diff > 0 && diff > fpepsilon)) {
+ G_message(_("EW extent does not match EW resolution: %g cells difference"),
+ diff);
+ }
+ if (cellhd->east - cellhd->west > 360.0) {
+ diff = (cellhd->east - cellhd->west - 360.0) / cellhd->ew_res;
+ if (diff > fpepsilon)
+ G_message(_("360 degree EW extent is exceeded by %g cells"),
+ diff);
+ }
+ else if (cellhd->east - cellhd->west < 360.0) {
+ diff = (360.0 - (cellhd->east - cellhd->west)) / cellhd->ew_res;
+ if (diff < 1.0 && diff > fpepsilon)
+ G_message(_("%g cells missing to cover 360 degrees"),
+ diff);
+ }
+
+ return lladjust;
+}
More information about the grass-commit
mailing list