v.distance: +geodesic (sync to trunk)

git-svn-id: https://svn.osgeo.org/grass/grass/branches/releasebranch_7_0@61966 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.distance/distance.c

@@ -56,9 +56,12 @@ int line2line(struct line_pnts *FPoints, int ftype,
     double tmp_dist, tmp_x, tmp_y, tmp_z, tmp_along;
     int ret = 1;
     static struct line_pnts *iPoints = NULL;
+    static struct line_pnts *LPoints = NULL;
     
-    if (!iPoints)
+    if (!iPoints) {
 	iPoints = Vect_new_line_struct();
+	LPoints = Vect_new_line_struct();
+    }
 
     *dist = PORT_DOUBLE_MAX;
 
@@ -75,12 +78,13 @@ int line2line(struct line_pnts *FPoints, int ftype,
     *ty = TPoints->y[0];
     *tz = TPoints->z[0];
 
+    tmp_z = 0;
+
     /* point -> point */
     if ((ftype & GV_POINTS) && (ttype & GV_POINTS)) {
 	Vect_line_distance(TPoints, FPoints->x[0], FPoints->y[0],
-			   FPoints->z[0], with_z, tx, ty,tz, dist, 
+			   FPoints->z[0], with_z, tx, ty, tz, dist, 
 			   NULL, talong);
-
     }
 
     /* point -> line and line -> line */
@@ -176,6 +180,45 @@ int line2line(struct line_pnts *FPoints, int ftype,
 	}
     }
 
+    if (geodesic) {
+	if (*fx != *tx || *fy != *ty || (with_z && *fz != *tz)) {
+	    Vect_reset_line(LPoints);
+	    Vect_append_point(LPoints, *fx, *fy, *fz);
+	    Vect_append_point(LPoints, *tx, *ty, *tz);
+	    *dist = Vect_line_geodesic_length(LPoints);
+	}
+	/* falong */
+	if (FPoints->x[0] != *fx || FPoints->y[0] != *fy ||
+	    (with_z && FPoints->z[0] != *fz)) {
+	
+	    fseg = Vect_line_distance(FPoints, *tx, *ty, *tz,
+			       with_z, &tmp_x, &tmp_y, &tmp_z,
+			       &tmp_dist, NULL, &tmp_along);
+
+	    Vect_reset_line(LPoints);
+	    for (i = 0; i < fseg; i++)
+		Vect_append_point(LPoints, FPoints->x[i], FPoints->y[i],
+		                  FPoints->z[i]);
+	    Vect_append_point(LPoints, *fx, *fy, *fz);
+	    *falong = Vect_line_geodesic_length(LPoints);
+	}
+	/* talong */
+	if (TPoints->x[0] != *tx || TPoints->y[0] != *ty ||
+	    (with_z && TPoints->z[0] != *tz)) {
+	
+	    tseg = Vect_line_distance(TPoints, *fx, *fy, *fz,
+			       with_z, &tmp_x, &tmp_y, &tmp_z,
+			       &tmp_dist, NULL, &tmp_along);
+
+	    Vect_reset_line(LPoints);
+	    for (i = 0; i < tseg; i++)
+		Vect_append_point(LPoints, TPoints->x[i], TPoints->y[i],
+		                  TPoints->z[i]);
+	    Vect_append_point(LPoints, *tx, *ty, *tz);
+	    *talong = Vect_line_geodesic_length(LPoints);
+	}
+    }
+
     return ret;
 }
 

vector/v.distance/main.c

@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
     struct Map_info From, To, Out, *Outp;
     int from_type, to_type, from_field, to_field, with_z;
     double max, min;
-    double *max_step;
+    double *max_step, max_map;
     int n_max_steps, curr_step;
     struct line_pnts *FPoints, *TPoints;
     struct line_cats *FCats, *TCats;
@@ -66,6 +66,7 @@ int main(int argc, char *argv[])
     double fx, fy, fz, falong, fangle, tx, ty, tz, talong, tangle, dist;
     double tmp_fx, tmp_fy, tmp_fz, tmp_falong, tmp_fangle, tmp_dist;
     double tmp_tx, tmp_ty, tmp_tz, tmp_talong, tmp_tangle;
+    int geodesic;
     struct field_info *Fi, *toFi;
     dbString stmt, dbstr;
     dbDriver *driver, *to_driver;
@@ -236,6 +237,8 @@ int main(int argc, char *argv[])
     if (flag.all->answer)
 	do_all = TRUE;
 
+    geodesic = G_projection() == PROJECTION_LL;
+
     /* Read upload and column options */
     /* count */
     i = 0;
@@ -295,7 +298,8 @@ int main(int argc, char *argv[])
 
     /* Open 'from' vector */
     Vect_set_open_level(2);
-    Vect_open_old2(&From, opt.from->answer, "", opt.from_field->answer);
+    if (Vect_open_old2(&From, opt.from->answer, "", opt.from_field->answer) < 0)
+	G_fatal_error(_("Unable to open vector map <%s>"), opt.from->answer);
 
     from_field = Vect_get_field_number(&From, opt.from_field->answer);
 
@@ -313,7 +317,8 @@ int main(int argc, char *argv[])
     
     /* Open 'to' vector */
     Vect_set_open_level(2);
-    Vect_open_old2(&To, opt.to->answer, "", opt.to_field->answer);
+    if (Vect_open_old2(&To, opt.to->answer, "", opt.to_field->answer) < 0)
+	G_fatal_error(_("Unable to open vector map <%s>"), opt.to->answer);
 
     ntolines = Vect_get_num_primitives(&To, to_type);
     ntoareas = 0;
@@ -333,7 +338,10 @@ int main(int argc, char *argv[])
 
     /* Open output vector */
     if (opt.out->answer) {
-	Vect_open_new(&Out, opt.out->answer, WITHOUT_Z);
+	if (Vect_open_new(&Out, opt.out->answer, WITHOUT_Z) < 0)
+	    G_fatal_error(_("Unable to create vector map <%s>"),
+			    opt.out->answer);
+
 	Vect_hist_command(&Out);
 	Outp = &Out;
     }
@@ -344,9 +352,10 @@ int main(int argc, char *argv[])
     /* TODO: add maxdist = -1 to Vect_select_ !!! */
     /* Calc maxdist */
     n_max_steps = 1;
+    max_map = max;
     if (max != 0) {
 	struct bound_box fbox, tbox;
-	double dx, dy, dz, tmp_max;
+	double dx, dy, dz;
 
 	Vect_get_map_box(&From, &fbox);
 	Vect_get_map_box(&To, &tbox);
@@ -355,21 +364,29 @@ int main(int argc, char *argv[])
 
 	dx = fbox.E - fbox.W;
 	dy = fbox.N - fbox.S;
+
 	if (Vect_is_3d(&From))
 	    dz = fbox.T - fbox.B;
 	else
 	    dz = 0.0;
 
-	tmp_max = sqrt(dx * dx + dy * dy + dz * dz);
-	if (max < 0)
-	    max = tmp_max;
+	max_map = sqrt(dx * dx + dy * dy + dz * dz);
+	if (max < 0) {
+	    if (geodesic)
+		max = PORT_DOUBLE_MAX;
+	    else
+		max = max_map;
+	}
 
 	/* how to determine a reasonable number of steps to increase the search box? */
 	/* with max > 0 but max <<< tmp_max, 2 steps are sufficient, first 0 then max
 	 * a reasonable number of steps also depends on the number of features in To
 	 * e.g. only one area in To, no need to step */
 
-	n_max_steps = sqrt(nto) * max / tmp_max;
+	if (geodesic)
+	    n_max_steps = sqrt(nto);
+	else
+	    n_max_steps = sqrt(nto) * max / max_map;
 	/* max 9 steps from testing */
 	if (n_max_steps > 9)
 	    n_max_steps = 9;
@@ -378,8 +395,8 @@ int main(int argc, char *argv[])
 	if (n_max_steps > nto)
 	    n_max_steps = nto;
 
-	G_debug(2, "max = %f", max);
-	G_debug(2, "maximum reasonable search distance = %f", tmp_max);
+	G_debug(2, "max = %g", max);
+	G_debug(2, "maximum reasonable search distance = %g", max_map);
 	G_debug(2, "n 'to' features = %d", nto);
 	G_debug(2, "n_max_steps = %d", n_max_steps);
     }
@@ -391,14 +408,17 @@ int main(int argc, char *argv[])
 	/* set up steps to increase search box */
 	max_step = G_malloc(n_max_steps * sizeof(double));
 	/* first step always 0 */
-	max_step[0] = (min < 0 ? 0 : min);
+	if (geodesic)
+	    max_step[0] = 0;
+	else
+	    max_step[0] = (min < 0 ? 0 : min);
 
 	for (curr_step = 1; curr_step < n_max_steps - 1; curr_step++) {
 	    /* for 9 steps, this would be max / [128, 64, 32, 16, 8, 4, 2] */
-	    max_step[curr_step] = max / (2 << (n_max_steps - 1 - curr_step));
+	    max_step[curr_step] = max_map / (2 << (n_max_steps - 1 - curr_step));
 	}
 	/* last step always max */
-	max_step[n_max_steps - 1] = max;
+	max_step[n_max_steps - 1] = max_map;
 	j = 1;
 	for (i = 1; i < n_max_steps; i++) {
 	    if (max_step[j - 1] < max_step[i]) {
@@ -410,7 +430,7 @@ int main(int argc, char *argv[])
     }
     else {
 	max_step = G_malloc(sizeof(double));
-	max_step[0] = max;
+	max_step[0] = max_map;
     }
 
     /* Open database driver */
@@ -608,6 +628,9 @@ int main(int argc, char *argv[])
 	    double box_edge = 0;
 	    int done = FALSE;
 
+	    if (geodesic)
+		tmp_min = 0;
+
 	    curr_step = 0;
 
 	    G_debug(3, "fline = %d", fline);
@@ -665,10 +688,10 @@ int main(int argc, char *argv[])
 		    }
 		}
 		else {
-		    box.E = fbox.E + max;
-		    box.W = fbox.W - max;
-		    box.N = fbox.N + max;
-		    box.S = fbox.S - max;
+		    box.E = fbox.E + max_map;
+		    box.W = fbox.W - max_map;
+		    box.N = fbox.N + max_map;
+		    box.S = fbox.S - max_map;
 		    box.T = PORT_DOUBLE_MAX;
 		    box.B = -PORT_DOUBLE_MAX;
 
@@ -684,7 +707,7 @@ int main(int argc, char *argv[])
 		    line2line(FPoints, ftype, TPoints, ttype,
 		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
-			      &tmp_dist, with_z, 0);
+			      &tmp_dist, with_z, geodesic);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
@@ -758,7 +781,7 @@ int main(int argc, char *argv[])
 		    line2area(&To, FPoints, ftype, aList->id[i], &aList->box[i],
 		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
-			      &tmp_dist, with_z, 0);
+			      &tmp_dist, with_z, geodesic);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
@@ -828,8 +851,17 @@ int main(int argc, char *argv[])
 		G_debug(4, "  dist = %f", dist);
 
 		if (!do_all && curr_step < n_max_steps) {
+		    double dist_map = dist;
+		    
+		    if (geodesic && tfeature > 0) {
+			double dx = fx - tx;
+			double dy = fy - ty;
+			double dz = fz - tz;
+
+			dist_map = sqrt(dx * dx + dy * dy + dz * dz);
+		    }
 		    /* enlarging the search box is possible */
-		    if (tfeature > 0 && dist > box_edge) {
+		    if (tfeature > 0 && dist_map > box_edge) {
 			/* line found but distance > search edge:
 			 * line bbox overlaps with search box, line itself is outside search box */
 			done = FALSE;
@@ -876,7 +908,10 @@ int main(int argc, char *argv[])
 	    double tmp_min = (min < 0 ? 0 : min);
 	    double box_edge = 0;
 	    int done = FALSE;
-	    
+
+	    if (geodesic)
+		tmp_min = 0;
+
 	    curr_step = 0;
 
 	    G_debug(3, "farea = %d", area);
@@ -935,10 +970,10 @@ int main(int argc, char *argv[])
 		    }
 		}
 		else {
-		    box.E = fbox.E + max;
-		    box.W = fbox.W - max;
-		    box.N = fbox.N + max;
-		    box.S = fbox.S - max;
+		    box.E = fbox.E + max_map;
+		    box.W = fbox.W - max_map;
+		    box.N = fbox.N + max_map;
+		    box.S = fbox.S - max_map;
 		    box.T = PORT_DOUBLE_MAX;
 		    box.B = -PORT_DOUBLE_MAX;
 
@@ -957,7 +992,7 @@ int main(int argc, char *argv[])
 		    line2area(&From, TPoints, ttype, area, &fbox,
 		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
 		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
-			      &tmp_dist, with_z, 0);
+			      &tmp_dist, with_z, geodesic);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
@@ -1043,7 +1078,7 @@ int main(int argc, char *argv[])
 		    poly = line2area(&From, TPoints, ttype, area, &fbox,
 		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
 		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
-			      &tmp_dist, with_z, 0);
+			      &tmp_dist, with_z, geodesic);
 
 		    if (poly == 3) {
 			/* 'to' outer ring is outside 'from' area,
@@ -1067,7 +1102,7 @@ int main(int argc, char *argv[])
 			    poly = line2area(&To, FPoints, ttype, tarea, &aList->box[i],
 				      &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 				      &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
-				      &tmp_dist, with_z, 0);
+				      &tmp_dist, with_z, geodesic);
 
 			    /* inside isle ? */
 			    poly = poly == 2;
@@ -1087,7 +1122,7 @@ int main(int argc, char *argv[])
 				line2area(&From, TPoints, ttype, area, &fbox,
 					  &tmp2_tx, &tmp2_ty, &tmp2_tz, &tmp2_talong, &tmp2_tangle,
 					  &tmp2_fx, &tmp2_fy, &tmp2_fz, &tmp2_falong, &tmp2_fangle,
-					  &tmp2_dist, with_z, 0);
+					  &tmp2_dist, with_z, geodesic);
 
 				if (tmp2_dist < tmp_dist) {
 				    tmp_dist = tmp2_dist;
@@ -1167,8 +1202,17 @@ int main(int argc, char *argv[])
 		}
 
 		if (!do_all && curr_step < n_max_steps) {
+		    double dist_map = dist;
+		    
+		    if (geodesic && tfeature > 0) {
+			double dx = fx - tx;
+			double dy = fy - ty;
+			double dz = fz - tz;
+
+			dist_map = sqrt(dx * dx + dy * dy + dz * dz);
+		    }
 		    /* enlarging the search box is possible */
-		    if (tfeature > 0 && dist > box_edge) {
+		    if (tfeature > 0 && dist_map > box_edge) {
 			/* area found but distance > search edge:
 			 * area bbox overlaps with search box, area itself is outside search box */
 			done = FALSE;
@@ -1262,7 +1306,10 @@ int main(int argc, char *argv[])
     else if (do_all && opt.table->answer) {	/* create new table */
 	db_set_string(&stmt, "create table ");
 	db_append_string(&stmt, opt.table->answer);
-	db_append_string(&stmt, " (from_cat integer");
+	if (Outp)
+	    db_append_string(&stmt, " (cat integer,from_cat integer");
+	else
+	    db_append_string(&stmt, " (from_cat integer");
 
 	j = 0;
 	while (Upload[j].upload != END) {
@@ -1331,6 +1378,8 @@ int main(int argc, char *argv[])
 	    Vect_reset_line(FPoints);
 	    Vect_reset_cats(FCats);
 
+	    Vect_cat_set(FCats, 1, i);
+
 	    Vect_append_point(FPoints, Near[i].from_x, Near[i].from_y, 0);
 
 	    if (Near[i].dist == 0) {
@@ -1385,8 +1434,13 @@ int main(int argc, char *argv[])
 	    }
 	}
 	else if (do_all) {		/* insert new record */
-	    sprintf(buf1, "insert into %s values ( %d ", opt.table->answer,
-		    Near[i].from_cat);
+	    if (!Outp)
+		sprintf(buf1, "insert into %s values ( %d ", opt.table->answer,
+			Near[i].from_cat);
+	    else
+		sprintf(buf1, "insert into %s values ( %d, %d ", opt.table->answer,
+			i, Near[i].from_cat);
+
 	    db_set_string(&stmt, buf1);
 
 	    j = 0;
@@ -1611,6 +1665,15 @@ int main(int argc, char *argv[])
 
     Vect_close(&From);
     if (Outp != NULL) {
+	if (do_all && opt.table->answer) {
+	    dbConnection connection;
+
+	    db_set_default_connection();
+	    db_get_connection(&connection);
+	    Vect_map_add_dblink(Outp, 1, NULL, opt.table->answer, "cat", 
+				connection.databaseName,
+				connection.driverName);
+	}
 	Vect_build(Outp);
 	Vect_close(Outp);
     }

vector/v.distance/v.distance.html

@@ -53,11 +53,9 @@ updated to NULL.
 <p>The upload <b>column</b>(s) must already exist. Create one with 
 <em><a href="v.db.addcolumn.html">v.db.addcolumn</a></em>.
 
-<!-- needs Vect_line_geodesic_distance()
 <p>In lat-long locations <em>v.distance</em> gives distances 
-(<em>dist</em> and <em>to_along</em>) not in degrees but in meters  
-calculated as geodesic distances on a sphere.
--->
+(<em>dist</em>, <em>from_along</em>, and <em>to_along</em>) not in 
+degrees but in meters calculated as geodesic distances on a sphere.
 
 <h2>EXAMPLES</h2>