coordinate transformation based on GCPs is done better by v.rectify

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@49954 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.transform/Makefile

@@ -6,7 +6,7 @@ PGM=v.transform
 include $(MODULE_TOPDIR)/include/Make/Module.make
 
 DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(GISDEP)
-LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB) $(TRANSLIB)
+LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
 EXTRA_INC = $(VECT_INC)
 EXTRA_CFLAGS = $(VECT_CFLAGS)
 

vector/v.transform/creat_trans.c

@@ -1,58 +0,0 @@
-/*
- ****************************************************************************
- *
- * MODULE:       v.transform
- * AUTHOR(S):    See other files as well...
- *               Eric G. Miller <egm2@jps.net>
- * PURPOSE:      To transform a vector layer's coordinates via a set of tie
- *               points.
- * COPYRIGHT:    (C) 2002 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *              License (>=v2). Read the file COPYING that comes with GRASS
- *              for details.
- *
- *****************************************************************************/
-
-/*
- *  create_transform_conversion () - main driver routine to prepare
- *    the transformation equation.
- *
- *  Written by the GRASS Team, 02/16/90, -mh.
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <grass/glocale.h>
-#include <grass/gis.h>
-#include <grass/vector.h>
-#include "trans.h"
-#include "local_proto.h"
-
-int create_transform_from_file(struct file_info *Coord)
-{
-    int status;
-    int n_points;
-
-    init_transform_arrays();
-
-    n_points = 0;
-    /*  Get the coordinates from the file.  */
-    if ((n_points = get_coor_from_file(Coord->fp)) < 0)
-	G_fatal_error(_("Error reading coordinates file"));
-
-    status = setup_transform(n_points);
-
-    if (status != ALL_OK) {
-	G_message(_("Number of points that have been entered [%d]"),
-		  n_points);
-	print_transform_error(status);
-	G_fatal_error(_("Error creating transformation"));
-    }
-
-    if (G_verbose() > G_verbose_std())
-	print_transform_resids(n_points);
-
-    return (0);
-
-}				/*  create_transform_from_file()  */

vector/v.transform/get_coor.c

@@ -1,46 +0,0 @@
-/*
- ****************************************************************************
- *
- * MODULE:       v.transform
- * AUTHOR(S):    See other files as well...
- *               Eric G. Miller <egm2@jps.net>
- * PURPOSE:      Read all the registration (map) coordinates in from the file
- * COPYRIGHT:    (C) 2002 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *              License (>=v2). Read the file COPYING that comes with GRASS
- *              for details.
- *
- *****************************************************************************/
-
-#include <stdio.h>
-#include "trans.h"
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-int get_coor_from_file(FILE * fp)
-{
-
-    int i;
-    char buff[128];
-
-    for (i = 0; i < MAX_COOR; i++) {
-
-	if (fgets(buff, sizeof(buff), fp) == NULL)
-	    break;
-
-	if (sscanf(buff, "%lf %lf %lf %lf", &ax[i], &ay[i], &bx[i], &by[i]) !=
-	    4) {
-	    /* comment or illegal line found */
-	    if (!buff[0] == '#')
-		G_fatal_error(_("Reading coordinates from file."));
-	    else
-		i--;		/* just comment found */
-	}
-	use[i] = 1;
-
-    }				/*  for (i)  */
-
-    return (i);
-
-}				/*    get_coor_from_file ()   */

vector/v.transform/local_proto.h

@@ -1,24 +1,4 @@
 
-/* command.c */
-int set_default_options(struct file_info *, struct file_info *,
-			struct file_info *, struct command_flags *);
-/* creat_trans.c */
-int create_transform_from_file(struct file_info *);
-
-/* get_coor.c */
-int get_coor_from_file(FILE *);
-
-/* main.c */
-int main(int, char *[]);
-
-/* print_trans.c */
-int print_transform_resids(int);
-
-/* setup_trans.c */
-int setup_transform(int);
-int init_transform_arrays(void);
-int print_transform_error(int);
-
 /* trans_digit.c */
-int transform_digit_file(struct Map_info *, struct Map_info *, int,
+int transform_digit_file(struct Map_info *, struct Map_info *,
 			 double, int, double *, char *, char **, int);

vector/v.transform/main.c

@@ -39,26 +39,13 @@
 #include "trans.h"
 #include "local_proto.h"
 
-double ax[MAX_COOR];	/*  current map   */
-double ay[MAX_COOR];
-
-double bx[MAX_COOR];	/*  map we are going to   */
-double by[MAX_COOR];
-
-int use[MAX_COOR];	/*  where the coordinate came from */
-double residuals[MAX_COOR];
-double rms;
-
-/*  this may be used in the future  */
-int reg_cnt;		/*  count of registered points */
-
 int main(int argc, char *argv[])
 {
-    struct file_info Current, Trans, Coord;
+    struct file_info Current, Trans;
 
     struct GModule *module;
 
-    struct Option *vold, *vnew, *pointsfile, *xshift, *yshift, *zshift,
+    struct Option *vold, *vnew, *xshift, *yshift, *zshift,
 	*xscale, *yscale, *zscale, *zrot, *columns, *table, *field;
     struct Flag *tozero_flag, *print_mat_flag, *swap_flag;
 
@@ -68,7 +55,7 @@ int main(int argc, char *argv[])
     struct bound_box box;
 
     double ztozero;
-    double trans_params[7];	// xshift, ..., xscale, ..., zrot
+    double trans_params[7];	/* xshift, ..., xscale, ..., zrot */
 
     /* columns */
     unsigned int i;
@@ -106,13 +93,6 @@ int main(int argc, char *argv[])
 
     vnew = G_define_standard_option(G_OPT_V_OUTPUT);
 
-    pointsfile = G_define_standard_option(G_OPT_F_INPUT);
-    pointsfile->key = "pointsfile";
-    pointsfile->required = NO;
-    pointsfile->label = _("ASCII file holding transform coordinates");
-    pointsfile->description = _("If not given, transformation parameters "
-				"(xshift, yshift, zshift, xscale, yscale, zscale, zrot) are used instead");
-
     xshift = G_define_option();
     xshift->key = "xshift";
     xshift->type = TYPE_DOUBLE;
@@ -207,20 +187,6 @@ int main(int argc, char *argv[])
 		       "Otherwise the table is overwritten."));
     }
 
-    if (pointsfile->answer != NULL) {
-	strcpy(Coord.name, pointsfile->answer);
-    }
-    else {
-	Coord.name[0] = '\0';
-    }
-
-    /* open coord file */
-    if (Coord.name[0] != '\0') {
-	if ((Coord.fp = fopen(Coord.name, "r")) == NULL)
-	    G_fatal_error(_("Unable to open file with coordinates <%s>"),
-			  Coord.name);
-    }
-
     /* open vector maps */
     Vect_open_old2(&Old, vold->answer, "", field->answer);
     Vect_open_new(&New, vnew->answer, Vect_is_3d(&Old) || tozero_flag->answer ||
@@ -247,14 +213,6 @@ int main(int argc, char *argv[])
     Vect_set_zone(&New, 0);
     Vect_set_thresh(&New, 0.0);
 
-    /* points file */
-    if (Coord.name[0]) {
-	create_transform_from_file(&Coord);
-
-	if (Coord.name[0] != '\0')
-	    fclose(Coord.fp);
-    }
-
     Vect_get_map_box(&Old, &box);
 
     /* z to zero */
@@ -312,7 +270,7 @@ int main(int argc, char *argv[])
     trans_params[IDX_ZROT] = atof(zrot->answer);
 
     G_important_message(_("Tranforming features..."));
-    transform_digit_file(&Old, &New, Coord.name[0] ? 1 : 0,
+    transform_digit_file(&Old, &New,
 			 ztozero, swap_flag->answer, trans_params,
 			 table->answer, columns_name, Vect_get_field_number(&Old, field->answer));
 
@@ -328,11 +286,7 @@ int main(int argc, char *argv[])
     G_verbose_message(_(" N: %-10.3f    S: %-10.3f"), box.N, box.S);
     G_verbose_message(_(" E: %-10.3f    W: %-10.3f"), box.E, box.W);
     G_verbose_message(_(" B: %6.3f    T: %6.3f"), box.B, box.T);
-    
-    /* print the transformation matrix if requested */
-    if (print_mat_flag->answer)
-      print_transform_matrix();
-    
+
     Vect_close(&New);
 
     G_done_msg(" ");

vector/v.transform/print_trans.c

@@ -1,49 +0,0 @@
-/*
- ****************************************************************************
- *
- * MODULE:       v.transform
- * AUTHOR(S):    See other files as well...
- *               Eric G. Miller <egm2@jps.net>
- * PURPOSE:      To transform a vector layer's coordinates via a set of tie
- *               points.
- * COPYRIGHT:    (C) 2002 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *              License (>=v2). Read the file COPYING that comes with GRASS
- *              for details.
- *
- *****************************************************************************/
-
-
-#include <stdio.h>
-#include <grass/gis.h>
-#include "trans.h"
-#include <grass/glocale.h>
-
-int print_transform_resids(int n_points)
-{
-    int i;
-
-    fprintf(stdout, "\n");
-    fprintf(stdout, "                CHECK MAP RESIDUALS\n");
-    fprintf(stdout, "                Current Map                  New Map\n");
-    fprintf(stdout,
-	    " POINT      X coord    Y coord  |        X coord   Y coord    |      residuals\n");
-    fprintf(stdout, "\n");
-
-    for (i = 0; i < n_points; i++) {
-
-	if (use[i])
-	    fprintf(stdout,
-		    " %2d.  %12.2f %12.2f | %12.2f   %12.2f | %12.2f\n",
-		    i + 1, ax[i], ay[i], bx[i], by[i], residuals[i]);
-
-    }
-
-    fprintf(stdout, "\n\n  Number of points: %d\n", n_points);
-    fprintf(stdout, "  Residual mean average: %f\n", rms);
-    fprintf(stdout, "\n");
-
-    return (0);
-
-}				/*  print_transform_resid()  */

vector/v.transform/setup_trans.c

@@ -1,94 +0,0 @@
-/*
- ****************************************************************************
- *
- * MODULE:       v.transform
- * AUTHOR(S):    See other files as well...
- *               Eric G. Miller <egm2@jps.net>
- * PURPOSE:      To transform a vector layer's coordinates via a set of tie
- *               points.
- * COPYRIGHT:    (C) 2002 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *              License (>=v2). Read the file COPYING that comes with GRASS
- *              for details.
- *
- *****************************************************************************/
-
-/*  
- *  Functions in this file:
- *     setup_transform ( n_points)
- *         returns:  ALL_OK, POINTS_NOT_SPREAD, NEED_MORE_POINTS
- *     init_transform_arrays()
- *     print_transform_error(stat)
- *
- *  Written by the GRASS Team, 02/16/90, -mh .
- */
-
-#include	<stdlib.h>
-#include	<stdio.h>
-#include <grass/glocale.h>
-#include <grass/gis.h>
-#include	"trans.h"
-
-int setup_transform(int n_points)
-{
-
-    int status;
-
-
-    /*  compute_transformation_coef() returns:
-     *   -2,  Not enough points
-     *   1,   everything is okay
-     *   -1,  points weren't spread out enough
-     */
-
-    /*  if there are enough points registered compute residuals  */
-    if (n_points >= MIN_COOR)
-	status = compute_transformation_coef(ax, ay, bx, by, use, MAX_COOR);
-    else
-	return (NEED_MORE_POINTS);
-
-    if (status != ALL_OK)
-	return (POINTS_NOT_SPREAD);
-
-    residuals_a_predicts_b(ax, ay, bx, by, use, MAX_COOR, residuals, &rms);
-
-    return (ALL_OK);
-
-
-}				/*  setup_transform()   */
-
-int init_transform_arrays(void)
-{
-    int i;
-
-    /*  initiliaze use[],  no points valid  */
-    for (i = 0; i < MAX_COOR; ++i) {
-	*(use + i) = 0;
-	*(bx + i) = 0;
-	*(by + i) = 0;
-	*(residuals + i) = 0;
-    }
-
-    reg_cnt = 0;
-
-    return 0;
-}
-
-
-int print_transform_error(int stat)
-{
-    switch (stat) {
-    case POINTS_NOT_SPREAD:
-	G_fatal_error(_("The points weren't spread out enough."));
-	break;
-    case NEED_MORE_POINTS:
-	G_fatal_error(_("You need to enter at least %d points."), MIN_COOR);
-	break;
-    default:
-	G_fatal_error("Your calling print_transform_error() with no error.");
-	break;
-    }
-
-    return 0;
-}

vector/v.transform/trans.h

@@ -12,15 +12,6 @@
 *  DEFINES:       *
 *******************/
 
-#define		MIN_COOR	4
-#define		MAX_COOR	1000
-
-
-/*  return status for setup_transform() and the transform library  */
-#define  ALL_OK  1
-#define  POINTS_NOT_SPREAD  -1
-#define  NEED_MORE_POINTS  -2
-
 #define  TRANS_MATRIX 0
 #define  TRANS_SHIFT  1
 
@@ -49,19 +40,6 @@
 * Yah, I made them global.  So shoot me.
 **/
 
-extern double ax[MAX_COOR];	/*  current map   */
-extern double ay[MAX_COOR];
-
-extern double bx[MAX_COOR];	/*  map we are going to   */
-extern double by[MAX_COOR];
-
-extern int use[MAX_COOR];	/*  where the coordinate came from */
-extern double residuals[MAX_COOR];
-extern double rms;
-
-/*  this may be used in the future  */
-extern int reg_cnt;		/*  count of registered points */
-
 
 /******************
 *  STRUCTS:       *
@@ -72,16 +50,5 @@ struct file_info
 {
     FILE *fp;
     char *mapset;
-    char name[80];
-    char full_name[256];
-};
-
-
-/* general flags that get set from the command line  */
-struct command_flags
-{
-    int verbose;		/*  do we print residual info  */
-    int usage;
+    char name[GPATH_MAX];
 };
-
-#include <grass/transform.h>

vector/v.transform/trans_digit.c

@@ -8,7 +8,7 @@
 #define PI M_PI
 
 int transform_digit_file(struct Map_info *Old, struct Map_info *New,
-			 int shift_file, double ztozero, int swap, double *trans_params_def,
+			 double ztozero, int swap, double *trans_params_def,
 			 char *table, char **columns, int field)
 {
     int i, type, cat, line;
@@ -124,28 +124,22 @@ int transform_digit_file(struct Map_info *Old, struct Map_info *New,
 
 	/* transform points */
 	for (i = 0; i < Points->n_points; i++) {
-	    if (shift_file) {
-		transform_a_into_b(Points->x[i], Points->y[i],
-				   &(Points->x[i]), &(Points->y[i]));
-	    }
-	    else {
-		G_debug(3, "idx=%d, cat=%d, xshift=%g, yshift=%g, zshift=%g, "
-			"xscale=%g, yscale=%g, zscale=%g, zrot=%g",
-			i, cat, trans_params[IDX_XSHIFT],
-			trans_params[IDX_YSHIFT], trans_params[IDX_ZSHIFT],
-			trans_params[IDX_XSCALE], trans_params[IDX_YSCALE],
-			trans_params[IDX_ZSCALE], trans_params[IDX_ZROT]);
-
-		/* transform point */
-		x = trans_params[IDX_XSHIFT] +
-		    trans_params[IDX_XSCALE] * Points->x[i] * cos(ang)
-		    - trans_params[IDX_YSCALE] * Points->y[i] * sin(ang);
-		y = trans_params[IDX_YSHIFT] +
-		    trans_params[IDX_XSCALE] * Points->x[i] * sin(ang)
-		    + trans_params[IDX_YSCALE] * Points->y[i] * cos(ang);
-		Points->x[i] = x;
-		Points->y[i] = y;
-	    }
+	    G_debug(3, "idx=%d, cat=%d, xshift=%g, yshift=%g, zshift=%g, "
+		    "xscale=%g, yscale=%g, zscale=%g, zrot=%g",
+		    i, cat, trans_params[IDX_XSHIFT],
+		    trans_params[IDX_YSHIFT], trans_params[IDX_ZSHIFT],
+		    trans_params[IDX_XSCALE], trans_params[IDX_YSCALE],
+		    trans_params[IDX_ZSCALE], trans_params[IDX_ZROT]);
+
+	    /* transform point */
+	    x = trans_params[IDX_XSHIFT] +
+		trans_params[IDX_XSCALE] * Points->x[i] * cos(ang)
+		- trans_params[IDX_YSCALE] * Points->y[i] * sin(ang);
+	    y = trans_params[IDX_YSHIFT] +
+		trans_params[IDX_XSCALE] * Points->x[i] * sin(ang)
+		+ trans_params[IDX_YSCALE] * Points->y[i] * cos(ang);
+	    Points->x[i] = x;
+	    Points->y[i] = y;
 
 	    /* ztozero shifts oldmap z to zero, zshift shifts rescaled object
 	     * to target elevation: */

vector/v.transform/v.transform.html

@@ -8,33 +8,16 @@ unreferenced vector maps or to modify existing geocoded maps.
 
 <h2>NOTES</h2>
 
-When using an ASCII table containing source and target coordinate pairs,
-in each row four coordinate values separated by white space have to be specified.
-Comments are permitted and have to be indicated by a '#' character.
-<p>Example for a points file of a linear transformation from XY to UTM coordinates
-(L: left, R: right, U: upper, L: lower, N, S, W, E):
-
-<div class="code"><pre>
-# Linear transformation from XY to UTM coordinates:
-# 4 maps corners defined
-# UL NW
-# UR NE
-# LR SW
-# LL SE
--584  585  598000 4920770
- 580  585  598020 4920770
- 580 -600  598020 4920750
--584 -600  598000 4920750
-</pre></div>
-
-<p>The ground control points may be also (ir)regularly distributed
-and can be more than four points.
+Coordinate transformation based on Ground Control Points (GCPs) is done 
+by <em><a href="v.rectify.html">v.rectify</a></em> and not supported by
+<em>v.transform</em>.
 
 <p>Transformation parameters (i.e. <em>xshift</em>, <em>yshift</em>,
 etc.) can be fetched from attribute table connected to the vector
 map. In this case vector objects can be transformed with different
 parameters based on their category number. If the parameter cannot be
 fetched from the table, default value is used instead.<p>
+
 <h3>Affine Transformation Matrix</h3>
 The affine transfomation matrix can optionally be printed with the '-m'
 flag. The format of the matrix is:
@@ -52,11 +35,12 @@ output of a similar operation performed in R.
 <h3>DXF/DWG drawings</h3>
 
 <p>Most DXF/DWG drawings are done within XY coordinate space. To transform 
-them to a national grid, we can use 'v.transform' with a 4 point 
-transformation.
+them to a national grid, we can use <em>v.transform</em> together with 
+<em>v.rectify</em> and a first-order 4 point transformation.
 
 <div class="code"><pre>
-v.transform -t in=watertowerXY out=watertowerUTM points=wt.points zscale=0.04 zshift=1320
+v.transform -t in=watertowerXY out=watertower_z zscale=0.04 zshift=1320
+v.rectify in=watertower_z out=watertowerUTM points=wt.points order=1
 </pre></div>
 
 <h3>Extrude 2D vector points to 3D based on attribute column values</h3>
@@ -103,6 +87,7 @@ The resulting map is a 3D vector map.
 
 <h2>SEE ALSO</h2>
 
+<em><a href="v.rectify.html">v.rectify</a></em>, 
 <em><a href="v.in.ogr.html">v.in.ogr</a></em>
 
 <h2>AUTHOR</h2>