[GRASS-SVN] r43886 - grass/trunk/imagery/i.rectify
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Oct 13 12:32:49 EDT 2010
Author: mmetz
Date: 2010-10-13 09:32:49 -0700 (Wed, 13 Oct 2010)
New Revision: 43886
Modified:
grass/trunk/imagery/i.rectify/get_wind.c
grass/trunk/imagery/i.rectify/global.h
grass/trunk/imagery/i.rectify/main.c
grass/trunk/imagery/i.rectify/rowcol.h
Log:
add resolution option, more robust resolution calculation
Modified: grass/trunk/imagery/i.rectify/get_wind.c
===================================================================
--- grass/trunk/imagery/i.rectify/get_wind.c 2010-10-13 13:23:05 UTC (rev 43885)
+++ grass/trunk/imagery/i.rectify/get_wind.c 2010-10-13 16:32:49 UTC (rev 43886)
@@ -1,3 +1,4 @@
+#include <math.h>
#include <grass/glocale.h>
#include "global.h"
#include "crs.h" /* CRS HEADER FILE */
@@ -2,11 +3,19 @@
-int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order)
+int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double res)
{
- double n, e;
+ double n, e, ad;
+ struct _corner {
+ double n, e;
+ } nw, ne, se, sw;
+ /* extends */
CRS_georef(w1->west, w1->north, &e, &n, E12, N12, order);
w2->north = w2->south = n;
w2->west = w2->east = e;
+ nw.n = n;
+ nw.e = e;
CRS_georef(w1->east, w1->north, &e, &n, E12, N12, order);
+ ne.n = n;
+ ne.e = e;
if (n > w2->north)
@@ -21,6 +30,8 @@
w2->west = e;
CRS_georef(w1->west, w1->south, &e, &n, E12, N12, order);
+ sw.n = n;
+ sw.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -31,6 +42,8 @@
w2->west = e;
CRS_georef(w1->east, w1->south, &e, &n, E12, N12, order);
+ se.n = n;
+ se.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -40,9 +53,63 @@
if (e < w2->west)
w2->west = e;
- w2->ns_res = (w2->north - w2->south) / w1->rows;
- w2->ew_res = (w2->east - w2->west) / w1->cols;
+ /* resolution */
+ if (res > 0)
+ w2->ew_res = w2->ns_res = res;
+ else {
+ /* this results in ugly res values, and ns_res != ew_res */
+ /* and is no good for rotation */
+ /*
+ w2->ns_res = (w2->north - w2->south) / w1->rows;
+ w2->ew_res = (w2->east - w2->west) / w1->cols;
+ */
+ /* alternative: account for rotation and order > 1 */
+
+ /* N-S extends along western and eastern edge */
+ w2->ns_res = (sqrt((nw.n - sw.n) * (nw.n - sw.n) +
+ (nw.e - sw.e) * (nw.e - sw.e)) +
+ sqrt((ne.n - se.n) * (ne.n - se.n) +
+ (ne.e - se.e) * (ne.e - se.e))) / (2.0 * w1->rows);
+
+ /* E-W extends along northern and southern edge */
+ w2->ew_res = (sqrt((nw.n - ne.n) * (nw.n - ne.n) +
+ (nw.e - ne.e) * (nw.e - ne.e)) +
+ sqrt((sw.n - se.n) * (sw.n - se.n) +
+ (sw.e - se.e) * (sw.e - se.e))) / (2.0 * w1->cols);
+
+ /* make ew_res = ns_res */
+ w2->ns_res = (w2->ns_res + w2->ew_res) / 2.0;
+ w2->ew_res = w2->ns_res;
+
+ /* nice round values */
+ if (w2->ns_res > 1) {
+ if (w2->ns_res < 10) {
+ /* round to first decimal */
+ w2->ns_res = (int)(w2->ns_res * 10 + 0.5) / 10.0;
+ w2->ew_res = w2->ns_res;
+ }
+ else {
+ /* round to integer */
+ w2->ns_res = (int)(w2->ns_res + 0.5);
+ w2->ew_res = w2->ns_res;
+ }
+ }
+ }
+
+ /* adjust extends */
+ ad = w2->north > 0 ? 0.5 : -0.5;
+ w2->north = (int) (ceil(w2->north / w2->ns_res) + ad) * w2->ns_res;
+ ad = w2->south > 0 ? 0.5 : -0.5;
+ w2->south = (int) (floor(w2->south / w2->ns_res) + ad) * w2->ns_res;
+ ad = w2->east > 0 ? 0.5 : -0.5;
+ w2->east = (int) (ceil(w2->east / w2->ew_res) + ad) * w2->ew_res;
+ ad = w2->west > 0 ? 0.5 : -0.5;
+ w2->west = (int) (floor(w2->west / w2->ew_res) + ad) * w2->ew_res;
+
+ w2->rows = (w2->north - w2->south + w2->ns_res / 2.0) / w2->ns_res;
+ w2->cols = (w2->east - w2->west + w2->ew_res / 2.0) / w2->ew_res;
+
G_message(_("Region N=%f S=%f E=%f W=%f"), w2->north, w2->south,
w2->east, w2->west);
G_message(_("Resolution EW=%f NS=%f"), w2->ew_res, w2->ns_res);
Modified: grass/trunk/imagery/i.rectify/global.h
===================================================================
--- grass/trunk/imagery/i.rectify/global.h 2010-10-13 13:23:05 UTC (rev 43885)
+++ grass/trunk/imagery/i.rectify/global.h 2010-10-13 16:32:49 UTC (rev 43886)
@@ -55,7 +55,7 @@
/* get_wind.c */
int get_target_window(int);
-int georef_window(struct Cell_head *, struct Cell_head *, int);
+int georef_window(struct Cell_head *, struct Cell_head *, int, double);
/* matrix.c */
int compute_georef_matrix(struct Cell_head *, struct Cell_head *, int);
Modified: grass/trunk/imagery/i.rectify/main.c
===================================================================
--- grass/trunk/imagery/i.rectify/main.c 2010-10-13 13:23:05 UTC (rev 43885)
+++ grass/trunk/imagery/i.rectify/main.c 2010-10-13 16:32:49 UTC (rev 43886)
@@ -69,10 +69,10 @@
char group[INAME_LEN], extension[INAME_LEN];
char result[NFILES][15];
int order; /* ADDED WITH CRS MODIFICATIONS */
- int n, i, m, k;
+ int n, i, m, k = 0;
int got_file = 0;
- struct Option *grp, *val, *ifile, *ext;
+ struct Option *grp, *val, *ifile, *ext, *tres;
struct Flag *c, *a;
struct GModule *module;
@@ -109,6 +109,12 @@
val->required = YES;
val->description = _("Rectification polynom order (1-3)");
+ tres = G_define_option();
+ tres->key = "res";
+ tres->type = TYPE_DOUBLE;
+ tres->required = NO;
+ tres->description = _("Target resolution (ignored if -c flag used)");
+
c = G_define_flag();
c->key = 'c';
c->description =
@@ -143,8 +149,11 @@
MAXORDER);
/* determine the number of files in this group */
- if (I_get_group_ref(group, &ref) <= 0)
+ if (I_get_group_ref(group, &ref) <= 0) {
+ G_warning(_("Location: %s"), G_location());
+ G_warning(_("Mapset: %s"), G_mapset());
G_fatal_error(_("Group <%s> does not exist"), grp->answer);
+ }
if (ref.nfiles <= 0) {
G_important_message(_("Group <%s> contains no raster maps; run i.group"),
@@ -184,21 +193,20 @@
/* get the target */
get_target(group);
-
- if (c->answer) {
- /* Use current Region */
- select_target_env();
- G__get_window(&target_window, "", "WIND", G_mapset());
- select_current_env();
- }
- else {
+ if (!c->answer) {
+ double res = -1;
+
+ if (tres->answer) {
+ if (!((res = atof(tres->answer)) > 0))
+ G_warning(_("Target resolution must be > 0, ignored"));
+ }
/* Calculate smallest region */
if (a->answer)
Rast_get_cellhd(ref.file[0].name, ref.file[0].mapset, &cellhd);
else
Rast_get_cellhd(ifile->answers[0], ref.file[0].mapset, &cellhd);
- georef_window(&cellhd, &target_window, order);
+ georef_window(&cellhd, &target_window, order, res);
}
G_verbose_message(_("Using region: N=%f S=%f, E=%f W=%f"), target_window.north,
Modified: grass/trunk/imagery/i.rectify/rowcol.h
===================================================================
--- grass/trunk/imagery/i.rectify/rowcol.h 2010-10-13 13:23:05 UTC (rev 43885)
+++ grass/trunk/imagery/i.rectify/rowcol.h 2010-10-13 16:32:49 UTC (rev 43886)
@@ -1 +1 @@
-#define ROWCOL short
+#define ROWCOL int
More information about the grass-commit
mailing list