[GRASS-SVN] r53677 - grass/trunk/raster/r.proj
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Nov 3 09:53:07 PDT 2012
Author: glynn
Date: 2012-11-03 09:53:06 -0700 (Sat, 03 Nov 2012)
New Revision: 53677
Modified:
grass/trunk/raster/r.proj/bordwalk.c
grass/trunk/raster/r.proj/main.c
grass/trunk/raster/r.proj/r.proj.h
Log:
Refactor bordwalk.c
Improve accuracy of -p/-g output
Modified: grass/trunk/raster/r.proj/bordwalk.c
===================================================================
--- grass/trunk/raster/r.proj/bordwalk.c 2012-11-03 16:50:04 UTC (rev 53676)
+++ grass/trunk/raster/r.proj/bordwalk.c 2012-11-03 16:53:06 UTC (rev 53677)
@@ -31,182 +31,203 @@
#include <grass/raster.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
+#include "r.proj.h"
-void bordwalk(struct Cell_head *from_hd, struct Cell_head *to_hd,
- struct pj_info *from_pj, struct pj_info *to_pj)
+#ifdef MIN
+#undef MIN
+#endif
+#define MIN(a, b) (((a) < (b)) ? (a) : (b))
+
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(a, b) (((a) > (b)) ? (a) : (b))
+
+static void debug(const char *name, const struct Cell_head *hd)
{
- double idx;
- double hx, hy;
- double xmin, xmax;
- double ymin, ymax;
+ G_debug(3, "%s: xmin: %f; xmax: %f; ymin: %f; ymax: %f",
+ name, hd->west, hd->east, hd->south, hd->north);
+}
- /* Set some (un)reasonable defaults before we walk the borders */
+static void update(struct Cell_head *to_hd, double hx, double hy)
+{
+ to_hd->east = MAX(to_hd->east, hx);
+ to_hd->west = MIN(to_hd->west, hx);
+ to_hd->north = MAX(to_hd->north, hy);
+ to_hd->south = MIN(to_hd->south, hy);
+}
- xmax = to_hd->west - 0.000001;
- xmin = to_hd->east + 0.000001;
- ymin = to_hd->north + 0.000001;
- ymax = to_hd->south - 0.000001;
+static void intersect(struct Cell_head *to_hd, const struct Cell_head *from_hd)
+{
+ to_hd->east = MIN(to_hd->east, from_hd->east);
+ to_hd->west = MAX(to_hd->west, from_hd->west);
+ to_hd->north = MIN(to_hd->north, from_hd->north);
+ to_hd->south = MAX(to_hd->south, from_hd->south);
+}
- /* Start walking */
+static int inside(const struct Cell_head *ref_hd, double hx, double hy)
+{
+ return
+ (hx <= ref_hd->east) &&
+ (hx >= ref_hd->west) &&
+ (hy <= ref_hd->north) &&
+ (hy >= ref_hd->south);
+}
+static void invert(struct Cell_head *cur_hd, const struct Cell_head *ref_hd,
+ double epsilon)
+{
+ cur_hd->east = ref_hd->west - epsilon;
+ cur_hd->west = ref_hd->east + epsilon;
+ cur_hd->north = ref_hd->south - epsilon;
+ cur_hd->south = ref_hd->north + epsilon;
+}
+
+static void proj_update(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ struct Cell_head *to_hd, double hx, double hy)
+{
+ if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
+ return;
+ update(to_hd, hx, hy);
+}
+
+void bordwalk1(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct Cell_head *from_hd, struct Cell_head *to_hd)
+{
+ double idx;
+
/* Top */
for (idx = from_hd->west + from_hd->ew_res / 2; idx < from_hd->east;
- idx += from_hd->ew_res) {
- hx = idx;
- hy = from_hd->north - from_hd->ns_res / 2;
- if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
- continue;
- /* check if we are within the region, but allow for some 'almost inside' points */
- /* (should probably be a factor based on input and output resolutions) */
- if (!(hx < to_hd->west - to_hd->ew_res) &&
- !(hx > to_hd->east + to_hd->ew_res) &&
- !(hy < to_hd->south - to_hd->ns_res) &&
- !(hy > to_hd->north + to_hd->ns_res)) {
- xmin = !(hx > xmin) ? hx : xmin;
- xmax = !(hx < xmax) ? hx : xmax;
- ymin = !(hy > ymin) ? hy : ymin;
- ymax = !(hy < ymax) ? hy : ymax;
- }
- }
+ idx += from_hd->ew_res)
+ proj_update(from_pj, to_pj, to_hd, idx, from_hd->north - from_hd->ns_res / 2);
- G_debug(3, "Top: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
- ymin, ymax);
+ debug("Top", to_hd);
/* Right */
for (idx = from_hd->north - from_hd->ns_res / 2; idx > from_hd->south;
- idx -= from_hd->ns_res) {
- hx = from_hd->east - from_hd->ew_res / 2;
- hy = idx;
- if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
- continue;
- if (!(hx < to_hd->west - to_hd->ew_res) &&
- !(hx > to_hd->east + to_hd->ew_res) &&
- !(hy < to_hd->south - to_hd->ns_res) &&
- !(hy > to_hd->north + to_hd->ns_res)) {
- xmin = !(hx > xmin) ? hx : xmin;
- xmax = !(hx < xmax) ? hx : xmax;
- ymin = !(hy > ymin) ? hy : ymin;
- ymax = !(hy < ymax) ? hy : ymax;
- }
- }
+ idx -= from_hd->ns_res)
+ proj_update(from_pj, to_pj, to_hd, from_hd->east - from_hd->ew_res / 2, idx);
- G_debug(3, "Right: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
- ymin, ymax);
+ debug("Right", to_hd);
/* Bottom */
for (idx = from_hd->east - from_hd->ew_res / 2; idx > from_hd->west;
- idx -= from_hd->ew_res) {
- hx = idx;
- hy = from_hd->south + from_hd->ns_res / 2;
- if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
- continue;
- if (!(hx < to_hd->west - to_hd->ew_res) &&
- !(hx > to_hd->east + to_hd->ew_res) &&
- !(hy < to_hd->south - to_hd->ns_res) &&
- !(hy > to_hd->north + to_hd->ns_res)) {
- xmin = !(hx > xmin) ? hx : xmin;
- xmax = !(hx < xmax) ? hx : xmax;
- ymin = !(hy > ymin) ? hy : ymin;
- ymax = !(hy < ymax) ? hy : ymax;
- }
- }
+ idx -= from_hd->ew_res)
+ proj_update(from_pj, to_pj, to_hd, idx, from_hd->south + from_hd->ns_res / 2);
- G_debug(3, "Bottom: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
- ymin, ymax);
+ debug("Bottom", to_hd);
/* Left */
for (idx = from_hd->south + from_hd->ns_res / 2; idx < from_hd->north;
- idx += from_hd->ns_res) {
- hx = from_hd->west + from_hd->ew_res / 2;
- hy = idx;
- if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
- continue;
- if (!(hx < to_hd->west - to_hd->ew_res) &&
- !(hx > to_hd->east + to_hd->ew_res) &&
- !(hy < to_hd->south - to_hd->ns_res) &&
- !(hy > to_hd->north + to_hd->ns_res)) {
- xmin = !(hx > xmin) ? hx : xmin;
- xmax = !(hx < xmax) ? hx : xmax;
- ymin = !(hy > ymin) ? hy : ymin;
- ymax = !(hy < ymax) ? hy : ymax;
- }
- }
+ idx += from_hd->ns_res)
+ proj_update(from_pj, to_pj, to_hd, from_hd->west + from_hd->ew_res / 2, idx);
- G_debug(3, "Left: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
- ymin, ymax);
+ debug("Left", to_hd);
+}
- /* check some special cases by reversing the projection */
+static int proj_inside(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct Cell_head *ref_hd, double hx, double hy)
+{
+ if (pj_do_proj(&hx, &hy, to_pj, from_pj) < 0)
+ return 0;
+ return inside(ref_hd, hx, hy);
+}
- if (xmin > to_hd->west) {
- hx = to_hd->west + to_hd->ew_res / 2;
- hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
- if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
- !(hx < from_hd->west) && !(hx > from_hd->east) &&
- !(hy < from_hd->south) && !(hy > from_hd->north))
- xmin = to_hd->west + to_hd->ew_res / 2;
+static void reverse_check(const struct pj_info *from_pj, const struct pj_info *to_pj,
+ const struct Cell_head *from_hd,
+ const struct Cell_head *to_hd,
+ struct Cell_head *cur_hd)
+{
+ if (cur_hd->west > to_hd->west) {
+ double hx = to_hd->west + to_hd->ew_res / 2;
+ double hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
+ if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ cur_hd->west = hx;
}
- if (xmax < to_hd->east) {
- hx = to_hd->east - to_hd->ew_res / 2;
- hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
- if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
- !(hx < from_hd->west) && !(hx > from_hd->east) &&
- !(hy < from_hd->south) && !(hy > from_hd->north))
- xmax = to_hd->east - to_hd->ew_res / 2;
+ if (cur_hd->east < to_hd->east) {
+ double hx = to_hd->east - to_hd->ew_res / 2;
+ double hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
+ if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ cur_hd->east = hx;
}
- if (ymin > to_hd->south) {
- hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
- hy = to_hd->south + to_hd->ns_res / 2;
- if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
- !(hx < from_hd->west) && !(hx > from_hd->east) &&
- !(hy < from_hd->south) && !(hy > from_hd->north))
- ymin = to_hd->south + to_hd->ns_res / 2;
+ if (cur_hd->south > to_hd->south) {
+ double hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
+ double hy = to_hd->south + to_hd->ns_res / 2;
+ if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ cur_hd->south = hy;
}
- if (ymax < to_hd->north) {
- hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
- hy = to_hd->north - to_hd->ns_res / 2;
- if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
- !(hx < from_hd->west) && !(hx > from_hd->east) &&
- !(hy < from_hd->south) && !(hy > from_hd->north))
- ymax = to_hd->north - to_hd->ns_res / 2;
+ if (cur_hd->north < to_hd->north) {
+ double hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
+ double hy = to_hd->north - to_hd->ns_res / 2;
+ if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+ cur_hd->north = hy;
}
+}
- G_debug(3, "Extra check: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin,
- xmax, ymin, ymax);
+static int outside(const struct Cell_head *cur_hd, const struct Cell_head *ref_hd)
+{
+ return
+ (cur_hd->west > ref_hd->east ) ||
+ (cur_hd->east < ref_hd->west ) ||
+ (cur_hd->south > ref_hd->north) ||
+ (cur_hd->north < ref_hd->south);
+}
+static void snap_to_grid(struct Cell_head *cur_hd, const struct Cell_head *ref_hd)
+{
+ int lidx = (int) floor(Rast_easting_to_col( cur_hd->west, ref_hd));
+ int ridx = (int) floor(Rast_easting_to_col( cur_hd->east, ref_hd));
+ int bidx = (int) floor(Rast_northing_to_row(cur_hd->south, ref_hd));
+ int tidx = (int) floor(Rast_northing_to_row(cur_hd->north, ref_hd));
+
+ cur_hd->west = Rast_col_to_easting( lidx + 0.0, ref_hd);
+ cur_hd->east = Rast_col_to_easting( ridx + 1.0, ref_hd);
+ cur_hd->south = Rast_row_to_northing(bidx + 1.0, ref_hd);
+ cur_hd->north = Rast_row_to_northing(tidx + 0.0, ref_hd);
+}
+
+void bordwalk(const struct Cell_head *from_hd, struct Cell_head *to_hd,
+ const struct pj_info *from_pj, const struct pj_info *to_pj)
+{
+ struct Cell_head cur_hd;
+
+ /* Set some (un)reasonable defaults before we walk the borders */
+
+ invert(&cur_hd, to_hd, 1.0e-6);
+
+ /* Start walking */
+
+ bordwalk1(from_pj, to_pj, from_hd, &cur_hd);
+
+ intersect(&cur_hd, to_hd);
+
+ /* check some special cases by reversing the projection */
+
+ reverse_check(from_pj, to_pj, from_hd, to_hd, &cur_hd);
+
+ debug("Extra check", &cur_hd);
+
/* if we still have some unresonable default minmax left, then abort */
- if ((xmin > to_hd->east) || (xmax < to_hd->west)
- || (ymin > to_hd->north) || (ymax < to_hd->south))
+ if (outside(&cur_hd, to_hd))
G_fatal_error(_("Input raster map is outside current region"));
- if (xmin < to_hd->west + to_hd->ew_res / 2)
- xmin = to_hd->west + to_hd->ew_res / 2;
- if (xmax > to_hd->east - to_hd->ew_res / 2)
- xmax = to_hd->east - to_hd->ew_res / 2;
- if (ymin < to_hd->south + to_hd->ns_res / 2)
- ymin = to_hd->south + to_hd->ns_res / 2;
- if (ymax > to_hd->north - to_hd->ns_res / 2)
- ymax = to_hd->north - to_hd->ns_res / 2;
+ intersect(&cur_hd, to_hd);
/* adjust to edges */
- idx = (int)floor(Rast_easting_to_col(xmin, to_hd));
- xmin = Rast_col_to_easting(idx + 0.0, to_hd);
- idx = (int)floor(Rast_easting_to_col(xmax, to_hd));
- xmax = Rast_col_to_easting(idx + 1.0, to_hd);
- idx = (int)floor(Rast_northing_to_row(ymin, to_hd));
- ymin = Rast_row_to_northing(idx + 1.0, to_hd);
- idx = (int)floor(Rast_northing_to_row(ymax, to_hd));
- ymax = Rast_row_to_northing(idx + 0.0, to_hd);
+ snap_to_grid(&cur_hd, to_hd);
- to_hd->west = (xmin < to_hd->west) ? to_hd->west : xmin;
- to_hd->east = (xmax > to_hd->east) ? to_hd->east : xmax;
- to_hd->south = (ymin < to_hd->south) ? to_hd->south : ymin;
- to_hd->north = (ymax > to_hd->north) ? to_hd->north : ymax;
+ intersect(to_hd, &cur_hd);
- G_debug(3, "Final check: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin,
- xmax, ymin, ymax);
+ debug("Final check", to_hd);
}
+
+void bordwalk2(const struct Cell_head *from_hd, struct Cell_head *to_hd,
+ const struct pj_info *from_pj, const struct pj_info *to_pj)
+{
+ bordwalk1(from_pj, to_pj, from_hd, to_hd);
+}
Modified: grass/trunk/raster/r.proj/main.c
===================================================================
--- grass/trunk/raster/r.proj/main.c 2012-11-03 16:50:04 UTC (rev 53676)
+++ grass/trunk/raster/r.proj/main.c 2012-11-03 16:53:06 UTC (rev 53677)
@@ -245,6 +245,7 @@
mapname = outmap->answer ? outmap->answer : inmap->answer;
if (mapname && !list->answer && !overwrite &&
+ !print_bounds->answer && !gprint_bounds->answer &&
G_find_raster(mapname, G_mapset()))
G_fatal_error(_("option <%s>: <%s> exists."), "output", mapname);
@@ -354,10 +355,15 @@
G_message(_("Input map <%s@%s> in location <%s>:"),
inmap->answer, setname, inlocation->answer);
- if (pj_do_proj(&iwest, &isouth, &iproj, &oproj) < 0)
- G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
- if (pj_do_proj(&ieast, &inorth, &iproj, &oproj) < 0)
- G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+ outcellhd.north = -1e9;
+ outcellhd.south = 1e9;
+ outcellhd.east = -1e9;
+ outcellhd.west = 1e9;
+ bordwalk2(&incellhd, &outcellhd, &iproj, &oproj);
+ inorth = outcellhd.north;
+ isouth = outcellhd.south;
+ ieast = outcellhd.east;
+ iwest = outcellhd.west;
G_format_northing(inorth, north_str, curr_proj);
G_format_northing(isouth, south_str, curr_proj);
Modified: grass/trunk/raster/r.proj/r.proj.h
===================================================================
--- grass/trunk/raster/r.proj/r.proj.h 2012-11-03 16:50:04 UTC (rev 53676)
+++ grass/trunk/raster/r.proj/r.proj.h 2012-11-03 16:53:06 UTC (rev 53677)
@@ -35,8 +35,10 @@
char *text; /* menu display - full description */
};
-extern void bordwalk(struct Cell_head *, struct Cell_head *, struct pj_info *,
- struct pj_info *);
+extern void bordwalk(const struct Cell_head *, struct Cell_head *,
+ const struct pj_info *, const struct pj_info *);
+extern void bordwalk2(const struct Cell_head *, struct Cell_head *,
+ const struct pj_info *, const struct pj_info *);
extern struct cache *readcell(int, const char *);
extern block *get_block(struct cache *, int);
extern void release_cache(struct cache *);
More information about the grass-commit
mailing list