[GRASS-SVN] r61966 - grass/branches/releasebranch_7_0/vector/v.distance
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Sep 14 23:30:32 PDT 2014
Author: mmetz
Date: 2014-09-14 23:30:32 -0700 (Sun, 14 Sep 2014)
New Revision: 61966
Modified:
grass/branches/releasebranch_7_0/vector/v.distance/distance.c
grass/branches/releasebranch_7_0/vector/v.distance/main.c
grass/branches/releasebranch_7_0/vector/v.distance/v.distance.html
Log:
v.distance: +geodesic (sync to trunk)
Modified: grass/branches/releasebranch_7_0/vector/v.distance/distance.c
===================================================================
--- grass/branches/releasebranch_7_0/vector/v.distance/distance.c 2014-09-15 05:59:33 UTC (rev 61965)
+++ grass/branches/releasebranch_7_0/vector/v.distance/distance.c 2014-09-15 06:30:32 UTC (rev 61966)
@@ -56,9 +56,12 @@
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 @@
*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 @@
}
}
+ 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;
}
Modified: grass/branches/releasebranch_7_0/vector/v.distance/main.c
===================================================================
--- grass/branches/releasebranch_7_0/vector/v.distance/main.c 2014-09-15 05:59:33 UTC (rev 61965)
+++ grass/branches/releasebranch_7_0/vector/v.distance/main.c 2014-09-15 06:30:32 UTC (rev 61966)
@@ -53,7 +53,7 @@
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 @@
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 @@
if (flag.all->answer)
do_all = TRUE;
+ geodesic = G_projection() == PROJECTION_LL;
+
/* Read upload and column options */
/* count */
i = 0;
@@ -295,7 +298,8 @@
/* 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 @@
/* 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 @@
/* 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 @@
/* 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 @@
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 @@
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 @@
/* 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 @@
}
else {
max_step = G_malloc(sizeof(double));
- max_step[0] = max;
+ max_step[0] = max_map;
}
/* Open database driver */
@@ -608,6 +628,9 @@
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 @@
}
}
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 @@
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 @@
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 @@
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 @@
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 @@
}
}
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 @@
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 @@
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 @@
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 @@
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 @@
}
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 @@
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 @@
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 @@
}
}
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 @@
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);
}
Modified: grass/branches/releasebranch_7_0/vector/v.distance/v.distance.html
===================================================================
--- grass/branches/releasebranch_7_0/vector/v.distance/v.distance.html 2014-09-15 05:59:33 UTC (rev 61965)
+++ grass/branches/releasebranch_7_0/vector/v.distance/v.distance.html 2014-09-15 06:30:32 UTC (rev 61966)
@@ -53,11 +53,9 @@
<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>
More information about the grass-commit
mailing list