v.distance: use Vect_line_geodesic_distance()

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

vector/v.distance/distance.c

@@ -1,9 +1,10 @@
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
 #include <grass/vector.h>
 #include "local_proto.h"
 
-/* TODO: geodesic distance for latlong */
+dist_func *line_distance;
 
 int get_line_box(const struct line_pnts *Points, 
                  struct bound_box *box)
@@ -37,6 +38,23 @@ int get_line_box(const struct line_pnts *Points,
     return 1;
 }
 
+/* segment angle */
+double sangle(struct line_pnts *Points, int segment)
+{
+    double dx, dy;
+
+    if (Points->n_points < 2 || segment < 1)
+	return -9;
+    if (segment >= Points->n_points)
+	G_fatal_error(_("Invalid segment number %d for %d points"),
+	              segment, Points->n_points);
+
+    dx = Points->x[segment] - Points->x[segment - 1];
+    dy = Points->y[segment] - Points->y[segment - 1];
+    
+    return atan2(dy, dx);
+}
+
 /* calculate distance parameters between two primitives
  * return 1 point to point
  * return 2 point to line
@@ -49,18 +67,15 @@ int line2line(struct line_pnts *FPoints, int ftype,
 	      double *tx, double *ty, double *tz,
 	      double *talong, double *tangle,
 	      double *dist,
-	      int with_z,
-	      int geodesic)
+	      int with_z)
 {
     int i, fseg, tseg, tmp_seg;
     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) {
 	iPoints = Vect_new_line_struct();
-	LPoints = Vect_new_line_struct();
     }
 
     *dist = PORT_DOUBLE_MAX;
@@ -82,7 +97,7 @@ int line2line(struct line_pnts *FPoints, int ftype,
 
     /* point -> point */
     if ((ftype & GV_POINTS) && (ttype & GV_POINTS)) {
-	Vect_line_distance(TPoints, FPoints->x[0], FPoints->y[0],
+	line_distance(TPoints, FPoints->x[0], FPoints->y[0],
 			   FPoints->z[0], with_z, tx, ty, tz, dist, 
 			   NULL, talong);
     }
@@ -90,11 +105,11 @@ int line2line(struct line_pnts *FPoints, int ftype,
     /* point -> line and line -> line */
     if ((ttype & GV_LINES)) {
 
-	tseg = 1;
+	fseg = tseg = 0;
 	/* calculate the min distance between each point in fline with tline */
 	for (i = 0; i < FPoints->n_points; i++) {
 
-	    tmp_seg = Vect_line_distance(TPoints, FPoints->x[i],
+	    tmp_seg = line_distance(TPoints, FPoints->x[i],
 	                                 FPoints->y[i], FPoints->z[i],
 					 with_z, &tmp_x, &tmp_y, &tmp_z,
 					 &tmp_dist, NULL, &tmp_along);
@@ -108,22 +123,39 @@ int line2line(struct line_pnts *FPoints, int ftype,
 		*tz = tmp_z;
 		*talong = tmp_along;
 		tseg = tmp_seg;
+		fseg = i + 1;
 	    }
 	}
-	Vect_point_on_line(TPoints, *talong, NULL, NULL, NULL,
-			   tangle, NULL);
+	*tangle = sangle(TPoints, tseg);
+
+	if (FPoints->n_points > 1 && fseg > 0) {
+	    int np = FPoints->n_points;
+	    
+	    fseg--;
+	    
+	    if (fseg > 0) {
+		FPoints->n_points = fseg + 1;
+		*falong = Vect_line_length(FPoints);
+		FPoints->n_points = np;
+	    }
+	    np = fseg;
+	    if (fseg == 0)
+		fseg = 1;
+	    *fangle = sangle(FPoints, fseg);
+	}
+
 	ret++;
     }
 
     /* line -> point and line -> line */
     if (ftype & GV_LINES) {
-	
-	fseg = 1;
+
+	fseg = tseg = 0;
 
 	/* calculate the min distance between each point in tline with fline */
 	for (i = 0; i < TPoints->n_points; i++) {
 
-	    tmp_seg = Vect_line_distance(FPoints, TPoints->x[i],
+	    tmp_seg = line_distance(FPoints, TPoints->x[i],
 			       TPoints->y[i], TPoints->z[i],
 			       with_z, &tmp_x, &tmp_y, &tmp_z,
 			       &tmp_dist, NULL, &tmp_along);
@@ -136,13 +168,28 @@ int line2line(struct line_pnts *FPoints, int ftype,
 		*tx = TPoints->x[i];
 		*ty = TPoints->y[i];
 		*tz = TPoints->z[i];
-		*talong = 0.;
-		*tangle = -9.;
 		fseg = tmp_seg;
+		tseg = i + 1;
+	    }
+	}
+	*fangle = sangle(FPoints, fseg);
+
+	if (TPoints->n_points > 1 && tseg > 0) {
+	    int np = TPoints->n_points;
+	    
+	    tseg--;
+
+	    if (tseg > 0) {
+		TPoints->n_points = tseg + 1;
+		*talong = Vect_line_length(TPoints);
+		TPoints->n_points = np;
 	    }
+	    np = tseg;
+	    if (tseg == 0)
+		tseg = 1;
+	    *tangle = sangle(TPoints, tseg);
 	}
-	Vect_point_on_line(FPoints, *falong, NULL, NULL, NULL,
-			   fangle, NULL);
+
 	ret++;
 	
 	if ((ttype & GV_LINES) && *dist > 0) {
@@ -162,63 +209,22 @@ int line2line(struct line_pnts *FPoints, int ftype,
 		    *fz = *tz = iPoints->z[0];
 		    
 		    /* falong, talong */
-		    Vect_line_distance(FPoints, iPoints->x[0],
+		    fseg = line_distance(FPoints, iPoints->x[0],
 				       iPoints->y[0], iPoints->z[0],
 				       with_z, NULL, NULL, NULL,
 				       NULL, NULL, falong);
-		    Vect_line_distance(TPoints, iPoints->x[0],
+		    tseg = line_distance(TPoints, iPoints->x[0],
 				       iPoints->y[0], iPoints->z[0],
 				       with_z, NULL, NULL, NULL,
 				       NULL, NULL, talong);
 		    /* fangle, tangle */
-		    Vect_point_on_line(FPoints, *falong, NULL, NULL, NULL,
-				       fangle, NULL);
-		    Vect_point_on_line(TPoints, *talong, NULL, NULL, NULL,
-				       tangle, NULL);
+		    *fangle = sangle(FPoints, fseg);
+		    *tangle = sangle(TPoints, tseg);
 		}
 	    }
 	}
     }
 
-    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;
 }
 
@@ -234,8 +240,7 @@ int line2area(const struct Map_info *To,
 	      double *tx, double *ty, double *tz,
 	      double *talong, double *tangle,
 	      double *dist,
-	      int with_z,
-	      int geodesic)
+	      int with_z)
 {
     int i, j;
     double tmp_dist;
@@ -306,7 +311,10 @@ int line2area(const struct Map_info *To,
 		line2line(Points, type, aPoints, GV_BOUNDARY,
 		          fx, fy, fz, falong, fangle,
 		          tx, ty, tz, talong, tangle,
-			  dist, with_z, geodesic);
+			  dist, with_z);
+
+		*talong = 0;
+		*tangle = -9;
 		
 		return 1;
 	    }
@@ -330,7 +338,7 @@ int line2area(const struct Map_info *To,
 			    line2line(Points, type, iPoints[j], GV_BOUNDARY,
 				      &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 				      &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
-				      &tmp_dist, with_z, geodesic);
+				      &tmp_dist, with_z);
 
 			    if (*dist > tmp_dist) {
 				*dist = tmp_dist;
@@ -344,7 +352,7 @@ int line2area(const struct Map_info *To,
 				*tx = tmp_tx;
 				*ty = tmp_ty;
 				*tz = tmp_tz;
-				*talong = tmp_talong;
+				*talong = 0;
 				*tangle = tmp_tangle;
 			    }
 
@@ -366,6 +374,9 @@ int line2area(const struct Map_info *To,
 		    *tx = Points->x[i];
 		    *ty = Points->y[i];
 		    *tz = Points->z[i];
+
+		    *fangle = *tangle = -9.;
+		    *falong = *talong = 0.;
 		    
 		    *dist = 0;
 
@@ -377,6 +388,10 @@ int line2area(const struct Map_info *To,
 		    if (*dist == 0) {
 			/* the line intersected with the isle boundary
 			 * -> line is partially inside the area */
+
+			*fangle = *tangle = -9.;
+			*falong = *talong = 0.;
+
 			return 1;
 		    }
 		    /* else continue with next point */
@@ -417,7 +432,9 @@ int line2area(const struct Map_info *To,
     line2line(Points, type, aPoints, GV_BOUNDARY,
 	      fx, fy, fz, falong, fangle,
 	      tx, ty, tz, talong, tangle,
-	      dist, with_z, geodesic);
+	      dist, with_z);
+
+    *talong = 0;
 
     if (*dist == 0)
 	return 1;

vector/v.distance/local_proto.h

@@ -1,6 +1,9 @@
 #include <grass/vector.h>
 #include <grass/dbmi.h>
 
+#ifndef _LOCAL_PROTO_
+#define _LOCAL_PROTO_
+
 /* define codes for characteristics of relation between two nearest features */
 #define CAT        1		/* category of nearest feature */
 #define FROM_X     2		/* x coordinate of nearest point on 'from' feature */
@@ -35,6 +38,12 @@ typedef struct
     char *column;		/* column name */
 } UPLOAD;
 
+
+typedef int dist_func(const struct line_pnts *, double, double, double, int,
+                       double *, double *, double *, double *, double *,
+                       double *);
+extern dist_func *line_distance;
+
 /* cmp.c */
 int cmp_near(const void *, const void *);
 int cmp_near_to(const void *, const void *);
@@ -50,8 +59,7 @@ int line2line(struct line_pnts *FPoints, int ftype,
 	      double *tx, double *ty, double *tz,
 	      double *talong, double *tangle,
 	      double *dist,
-	      int with_z,
-	      int geodesic);
+	      int with_z);
 int line2area(const struct Map_info *To,
 	      struct line_pnts *Points, int type,
 	      int area, const struct bound_box *abox,
@@ -60,8 +68,9 @@ int line2area(const struct Map_info *To,
 	      double *tx, double *ty, double *tz,
 	      double *talong, double *tangle,
 	      double *dist,
-	      int with_z,
-	      int geodesic);
+	      int with_z);
 
 /* print.c */
 int print_upload(NEAR *, UPLOAD *, int, dbCatValArray *, dbCatVal *, char *);
+
+#endif

vector/v.distance/main.c

@@ -238,6 +238,10 @@ int main(int argc, char *argv[])
 	do_all = TRUE;
 
     geodesic = G_projection() == PROJECTION_LL;
+    if (geodesic)
+	line_distance = Vect_line_geodesic_distance;
+    else
+	line_distance = Vect_line_distance;
 
     /* Read upload and column options */
     /* count */
@@ -707,7 +711,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, geodesic);
+			      &tmp_dist, with_z);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
@@ -781,7 +785,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, geodesic);
+			      &tmp_dist, with_z);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
@@ -992,7 +996,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, geodesic);
+			      &tmp_dist, with_z);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
@@ -1078,7 +1082,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, geodesic);
+			      &tmp_dist, with_z);
 
 		    if (poly == 3) {
 			/* 'to' outer ring is outside 'from' area,
@@ -1102,7 +1106,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, geodesic);
+				      &tmp_dist, with_z);
 
 			    /* inside isle ? */
 			    poly = poly == 2;
@@ -1122,7 +1126,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, geodesic);
+					  &tmp2_dist, with_z);
 
 				if (tmp2_dist < tmp_dist) {
 				    tmp_dist = tmp2_dist;
@@ -1306,7 +1310,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) {
@@ -1375,6 +1382,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) {
@@ -1429,8 +1438,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;
@@ -1655,6 +1669,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);
     }