[GRASS-SVN] r70656 - in grass/trunk: include/defs lib/gis
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Feb 22 01:02:27 PST 2017
Author: mmetz
Date: 2017-02-22 01:02:27 -0800 (Wed, 22 Feb 2017)
New Revision: 70656
Modified:
grass/trunk/include/defs/gis.h
grass/trunk/lib/gis/adj_cellhd.c
Log:
libgis: new G_adjust_window_ll() to automatically fix lat/lon errors
Modified: grass/trunk/include/defs/gis.h
===================================================================
--- grass/trunk/include/defs/gis.h 2017-02-22 01:42:33 UTC (rev 70655)
+++ grass/trunk/include/defs/gis.h 2017-02-22 09:02:27 UTC (rev 70656)
@@ -75,6 +75,7 @@
/* adj_cellhd.c */
void G_adjust_Cell_head(struct Cell_head *, int, int);
void G_adjust_Cell_head3(struct Cell_head *, int, int, int);
+int G_adjust_window_ll(struct Cell_head *cellhd);
/* alloc.c */
#define G_incr_void_ptr(ptr, size) \
Modified: grass/trunk/lib/gis/adj_cellhd.c
===================================================================
--- grass/trunk/lib/gis/adj_cellhd.c 2017-02-22 01:42:33 UTC (rev 70655)
+++ grass/trunk/lib/gis/adj_cellhd.c 2017-02-22 09:02:27 UTC (rev 70656)
@@ -16,6 +16,8 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+#define LL_TOLERANCE 10
+
/* TODO: find good thresholds */
/* deviation measured in cells */
static double llepsilon = 0.01;
@@ -103,16 +105,16 @@
/* (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"));
+ if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
+ G_verbose_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 (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
+ G_verbose_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"));
+ if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
+ G_verbose_message(_("NS and EW resolutions are different"));
ll_check_ns(cellhd);
ll_check_ew(cellhd);
@@ -243,16 +245,16 @@
/* (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"));
+ if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
+ G_verbose_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 (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
+ G_verbose_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"));
+ if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
+ G_verbose_message(_("NS and EW resolutions are different"));
ll_check_ns(cellhd);
ll_check_ew(cellhd);
@@ -305,18 +307,18 @@
}
/* very liberal thresholds */
- if (cellhd->north > 100.0)
+ if (cellhd->north > 90.0 + LL_TOLERANCE)
G_fatal_error(_("Illegal latitude for North: %g"), cellhd->north);
- if (cellhd->south < -100.0)
+ if (cellhd->south < -90.0 - LL_TOLERANCE)
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) {
+ if (cellhd->west < -360.0 - LL_TOLERANCE) {
G_debug(1, "East: %g", cellhd->east);
G_fatal_error(_("Illegal longitude for West: %g"), cellhd->west);
}
- if (cellhd->east > 370.0) {
+ if (cellhd->east > 360.0 + LL_TOLERANCE) {
G_debug(1, "West: %g", cellhd->west);
G_fatal_error(_("Illegal longitude for East: %g"), cellhd->east);
}
@@ -346,7 +348,7 @@
diff -= ncells;
if ((diff < 0 && diff < -fpepsilon) ||
(diff > 0 && diff > fpepsilon)) {
- G_message(_("NS extent does not match NS resolution: %g cells difference"),
+ G_verbose_message(_("NS extent does not match NS resolution: %g cells difference"),
diff);
}
@@ -355,10 +357,10 @@
if (diff < 0)
diff = -diff;
if (cellhd->north < 90.0 && diff < 1.0 ) {
- G_message(_("%g cells missing to reach 90 degree north"),
+ G_verbose_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)"),
+ G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
cellhd->north - 90.0);
/* check only, do not modify
cellhd->north = 90.0;
@@ -368,11 +370,11 @@
}
if (cellhd->north > 90.0) {
if (diff <= 0.5 + llepsilon) {
- G_message(_("90 degree north is exceeded by %g cells"),
+ G_verbose_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)"),
+ G_verbose_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);
@@ -386,7 +388,7 @@
if (diff < 0)
diff = -diff;
if (diff < llepsilon && diff > fpepsilon) {
- G_message(_("Subtle input data rounding error of north boundary (%g)"),
+ G_verbose_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);
@@ -405,10 +407,10 @@
if (diff < 0)
diff = -diff;
if (cellhd->south > -90.0 && diff < 1.0 ) {
- G_message(_("%g cells missing to reach 90 degree south"),
+ G_verbose_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)"),
+ G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
cellhd->south + 90.0);
/* check only, do not modify
cellhd->south = -90.0;
@@ -418,11 +420,11 @@
}
if (cellhd->south < -90.0) {
if (diff <= 0.5 + llepsilon) {
- G_message(_("90 degree south is exceeded by %g cells"),
+ G_verbose_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)"),
+ G_verbose_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);
@@ -436,7 +438,7 @@
if (diff < 0)
diff = -diff;
if (diff < llepsilon && diff > fpepsilon) {
- G_message(_("Subtle input data rounding error of south boundary (%g)"),
+ G_verbose_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);
@@ -476,21 +478,374 @@
diff -= ncells;
if ((diff < 0 && diff < -fpepsilon) ||
(diff > 0 && diff > fpepsilon)) {
- G_message(_("EW extent does not match EW resolution: %g cells difference"),
+ G_verbose_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"),
+ G_verbose_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"),
+ G_verbose_message(_("%g cells missing to cover 360 degree EW extent"),
diff);
}
return lladjust;
}
+
+/*!
+ * \brief Adjust window for lat/lon.
+ *
+ * This function tries to automatically fix fp precision issues and
+ * adjust rounding errors for lat/lon.
+ *
+ * <b>Note:</b> 3D values are not adjusted.
+ *
+ * \param[in,out] cellhd pointer to Cell_head structure
+ * \return 1 if window was adjusted
+ * \return 0 if window was not adjusted
+ */
+int G_adjust_window_ll(struct Cell_head *cellhd)
+{
+ int ll_adjust, res_adj;
+ double dsec, dsec2;
+ char buf[100], buf2[100];
+ double diff;
+ double old, new;
+ struct Cell_head cellhds; /* everyting in seconds, not degrees */
+
+ /* lat/lon checks */
+ if (cellhd->proj != PROJECTION_LL)
+ return 0;
+
+ /* put everything through ll_format + ll_scan */
+ G_llres_format(cellhd->ns_res, buf);
+ if (G_llres_scan(buf, &new) != 1)
+ G_fatal_error(_("Invalid NS resolution"));
+ cellhd->ns_res = new;
+
+ G_llres_format(cellhd->ew_res, buf);
+ if (G_llres_scan(buf, &new) != 1)
+ G_fatal_error(_("Invalid EW resolution"));
+ cellhd->ew_res = new;
+
+ G_lat_format(cellhd->north, buf);
+ if (G_lat_scan(buf, &new) != 1)
+ G_fatal_error(_("Invalid North"));
+ cellhd->north = new;
+
+ G_lat_format(cellhd->south, buf);
+ if (G_lat_scan(buf, &new) != 1)
+ G_fatal_error(_("Invalid South"));
+ cellhd->south = new;
+
+ G_lon_format(cellhd->west, buf);
+ if (G_lon_scan(buf, &new) != 1)
+ G_fatal_error(_("Invalid West"));
+ cellhd->west = new;
+
+ G_lon_format(cellhd->east, buf);
+ if (G_lon_scan(buf, &new) != 1)
+ G_fatal_error(_("Invalid East"));
+ cellhd->east = new;
+
+ /* convert to seconds */
+ cellhds = *cellhd;
+
+ old = cellhds.ns_res * 3600;
+ sprintf(buf, "%f", old);
+ sscanf(buf, "%lf", &new);
+ cellhds.ns_res = new;
+
+ old = cellhds.ew_res * 3600;
+ sprintf(buf, "%f", old);
+ sscanf(buf, "%lf", &new);
+ cellhds.ew_res = new;
+
+ old = cellhds.north * 3600;
+ sprintf(buf, "%f", old);
+ sscanf(buf, "%lf", &new);
+ cellhds.north = new;
+
+ old = cellhds.south * 3600;
+ sprintf(buf, "%f", old);
+ sscanf(buf, "%lf", &new);
+ cellhds.south = new;
+
+ old = cellhds.west * 3600;
+ sprintf(buf, "%f", old);
+ sscanf(buf, "%lf", &new);
+ cellhds.west = new;
+
+ old = cellhds.east * 3600;
+ sprintf(buf, "%f", old);
+ sscanf(buf, "%lf", &new);
+ cellhds.east = new;
+
+ ll_adjust = 0;
+
+ /* N - S */
+ /* resolution */
+ res_adj = 0;
+ old = cellhds.ns_res;
+
+ if (old > 0.4) {
+ /* round to nearest 0.1 sec */
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = fabs(dsec2 - dsec) / dsec;
+ if (diff > 0 && diff < llepsilon) {
+ G_llres_format(old / 3600, buf);
+ G_llres_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("NS resolution rounded from %s to %s"),
+ buf, buf2);
+ ll_adjust = 1;
+ res_adj = 1;
+ cellhds.ns_res = new;
+ }
+ }
+
+ if (res_adj) {
+ double n_off, s_off;
+
+ old = cellhds.north;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
+ n_off = diff;
+
+ old = cellhds.south;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
+ s_off = diff;
+
+ if (n_off < llepsilon || n_off <= s_off) {
+ old = cellhds.north;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = n_off;
+ if (diff > 0 && diff < llepsilon) {
+ G_lat_format(old / 3600, buf);
+ G_lat_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("North rounded from %s to %s"),
+ buf, buf2);
+ cellhds.north = new;
+ }
+
+ old = cellhds.south;
+ new = cellhds.north - cellhds.ns_res * cellhds.rows;
+ diff = fabs(new - old) / cellhds.ns_res;
+ if (diff > 0) {
+ G_lat_format(old / 3600, buf);
+ G_lat_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("South adjusted from %s to %s"),
+ buf, buf2);
+ }
+ cellhds.south = new;
+ }
+ else {
+ old = cellhds.south;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = s_off;
+ if (diff > 0 && diff < llepsilon) {
+ G_lat_format(old / 3600, buf);
+ G_lat_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("South rounded from %s to %s"),
+ buf, buf2);
+ cellhds.south = new;
+ }
+
+ old = cellhds.north;
+ new = cellhds.south + cellhds.ns_res * cellhds.rows;
+ diff = fabs(new - old) / cellhds.ns_res;
+ if (diff > 0) {
+ G_lat_format(old / 3600, buf);
+ G_lat_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("North adjusted from %s to %s"),
+ buf, buf2);
+ }
+ cellhds.north = new;
+ }
+ }
+ else {
+ old = cellhds.north;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
+ if (diff > 0 && diff < llepsilon) {
+ G_lat_format(old / 3600, buf);
+ G_lat_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("North rounded from %s to %s"),
+ buf, buf2);
+ ll_adjust = 1;
+ cellhds.north = new;
+ }
+
+ old = cellhds.south;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
+ if (diff > 0 && diff < llepsilon) {
+ G_lat_format(old / 3600, buf);
+ G_lat_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("South rounded from %s to %s"),
+ buf, buf2);
+ ll_adjust = 1;
+ cellhds.south = new;
+ }
+ }
+ cellhds.ns_res = (cellhds.north - cellhds.south) / cellhds.rows;
+
+ /* E - W */
+ /* resolution */
+ res_adj = 0;
+ old = cellhds.ew_res;
+
+ if (old > 0.4) {
+ /* round to nearest 0.1 sec */
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = fabs(dsec2 - dsec) / dsec;
+ if (diff > 0 && diff < llepsilon) {
+ G_llres_format(old / 3600, buf);
+ G_llres_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("EW resolution rounded from %s to %s"),
+ buf, buf2);
+ ll_adjust = 1;
+ res_adj = 1;
+ cellhds.ew_res = new;
+ }
+ }
+
+ if (res_adj) {
+ double w_off, e_off;
+
+ old = cellhds.west;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
+ w_off = diff;
+
+ old = cellhds.east;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
+ e_off = diff;
+
+ if (w_off < llepsilon || w_off <= e_off) {
+ old = cellhds.west;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = w_off;
+ if (diff > 0 && diff < llepsilon) {
+ G_lon_format(old / 3600, buf);
+ G_lon_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("West rounded from %s to %s"),
+ buf, buf2);
+ cellhds.west = new;
+ }
+
+ old = cellhds.east;
+ new = cellhds.west + cellhds.ew_res * cellhds.cols;
+ diff = fabs(new - old) / cellhds.ew_res;
+ if (diff > 0) {
+ G_lon_format(old / 3600, buf);
+ G_lon_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("East adjusted from %s to %s"),
+ buf, buf2);
+ }
+ cellhds.east = new;
+ }
+ else {
+ old = cellhds.east;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = e_off;
+ if (diff > 0 && diff < llepsilon) {
+ G_lon_format(old / 3600, buf);
+ G_lon_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("East rounded from %s to %s"),
+ buf, buf2);
+ cellhds.east = new;
+ }
+
+ old = cellhds.west;
+ new = cellhds.east - cellhds.ew_res * cellhds.cols;
+ diff = fabs(new - cellhds.west) / cellhds.ew_res;
+ if (diff > 0) {
+ G_lon_format(old / 3600, buf);
+ G_lon_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("West adjusted from %s to %s"),
+ buf, buf2);
+ }
+ cellhds.west = new;
+ }
+ }
+ else {
+ old = cellhds.west;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
+ if (diff > 0 && diff < llepsilon) {
+ G_lon_format(old / 3600, buf);
+ G_lon_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("West rounded from %s to %s"),
+ buf, buf2);
+ ll_adjust = 1;
+ cellhds.west = new;
+ }
+
+ old = cellhds.east;
+ dsec = old * 10;
+ dsec2 = floor(dsec + 0.5);
+ new = dsec2 / 10;
+ diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
+ if (diff > 0 && diff < llepsilon) {
+ G_lon_format(old / 3600, buf);
+ G_lon_format(new / 3600, buf2);
+ if (strcmp(buf, buf2))
+ G_verbose_message(_("East rounded from %s to %s"),
+ buf, buf2);
+ ll_adjust = 1;
+ cellhds.east = new;
+ }
+ }
+ cellhds.ew_res = (cellhds.east - cellhds.west) / cellhds.cols;
+
+ cellhd->ns_res = cellhds.ns_res / 3600;
+ cellhd->ew_res = cellhds.ew_res / 3600;
+ cellhd->north = cellhds.north / 3600;
+ cellhd->south = cellhds.south / 3600;
+ cellhd->west = cellhds.west / 3600;
+ cellhd->east = cellhds.east / 3600;
+
+ return ll_adjust;
+}
More information about the grass-commit
mailing list