[GRASS-SVN] r53055 - grass/branches/develbranch_6/raster/r.grow.distance
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Sep 2 05:52:54 PDT 2012
Author: neteler
Date: 2012-09-02 05:52:53 -0700 (Sun, 02 Sep 2012)
New Revision: 53055
Modified:
grass/branches/develbranch_6/raster/r.grow.distance/description.html
grass/branches/develbranch_6/raster/r.grow.distance/main.c
Log:
r.grow.distance: backport of geodetic distance support (trunk, r34429); +example
Modified: grass/branches/develbranch_6/raster/r.grow.distance/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.grow.distance/description.html 2012-09-02 12:49:21 UTC (rev 53054)
+++ grass/branches/develbranch_6/raster/r.grow.distance/description.html 2012-09-02 12:52:53 UTC (rev 53055)
@@ -4,7 +4,6 @@
distance to the nearest non-null cell in the input map and/or the
value of the nearest non-null cell.
-
<h2>NOTES</h2>
The user has the option of specifying four different metrics which
control the geometry in which grown cells are created, (controlled by
@@ -12,13 +11,14 @@
<i>Manhattan</i>, and <i>Maximum</i>.
<p>
-
The <i>Euclidean distance</i> or <i>Euclidean metric</i> is the "ordinary" distance
between two points that one would measure with a ruler, which can be
proven by repeated application of the Pythagorean theorem.
The formula is given by:
-<div class="code"><pre>d(dx,dy) = sqrt(dx^2 + dy^2)</pre></div>
+<div class="code"><pre>
+d(dx,dy) = sqrt(dx^2 + dy^2)
+</pre></div>
Cells grown using this metric would form isolines of distance that are
circular from a given point, with the distance given by the <b>radius</b>.
@@ -29,7 +29,6 @@
and is sufficient if only relative values are required.
<p>
-
The <i>Manhattan metric</i>, or <i>Taxicab geometry</i>, is a form of geometry in
which the usual metric of Euclidean geometry is replaced by a new
metric in which the distance between two points is the sum of the (absolute)
@@ -39,39 +38,48 @@
points' distance in taxicab geometry.
The formula is given by:
-<div class="code"><pre>d(dx,dy) = abs(dx) + abs(dy)</pre></div>
+<div class="code"><pre>
+d(dx,dy) = abs(dx) + abs(dy)
+</pre></div>
where cells grown using this metric would form isolines of distance that are
rhombus-shaped from a given point.
<p>
-
The <i>Maximum metric</i> is given by the formula
-<div class="code"><pre>d(dx,dy) = max(abs(dx),abs(dy))</pre></div>
+<div class="code"><pre>
+d(dx,dy) = max(abs(dx),abs(dy))
+</pre></div>
where the isolines of distance from a point are squares.
<h2>EXAMPLE</h2>
-Spearfish sample dataset
+Distance from the streams network (North Carolina sample dataset):
<div class="code"><pre>
-r.grow.distance input=roads distance=dist_from_roads
+g.region rast=streams_derived -p
+r.grow.distance input=streams_derived distance=dist_from_streams
</pre></div>
+<p>
+Distance from sea in meters in latitude-longitude location:
+<div class="code"><pre>
+g.region rast=sea -p
+r.grow.distance -m input=sea distance=dist_from_sea_geodetic metric=geodesic
+</pre></div>
<h2>SEE ALSO</h2>
<em>
-<a href="r.grow.html">r.grow</a><br>
-<a href="r.buffer.html">r.buffer</a><br>
-<a href="r.cost.html">r.cost</a><br>
+<a href="r.grow.html">r.grow</a>,
+<a href="r.buffer.html">r.buffer</a>,
+<a href="r.cost.html">r.cost</a>,
<a href="r.patch.html">r.patch</a>
</em>
<p>
-
<em>
<a href="http://en.wikipedia.org/wiki/Euclidean_metric">Wikipedia Entry:
Euclidean Metric</a><br>
Modified: grass/branches/develbranch_6/raster/r.grow.distance/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.grow.distance/main.c 2012-09-02 12:49:21 UTC (rev 53054)
+++ grass/branches/develbranch_6/raster/r.grow.distance/main.c 2012-09-02 12:52:53 UTC (rev 53055)
@@ -1,7 +1,7 @@
/****************************************************************************
*
- * MODULE: r.grow2
+ * MODULE: r.grow.distance
*
* AUTHOR(S): Marjorie Larson - CERL
* Glynn Clements
@@ -26,6 +26,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+static struct Cell_head window;
static int nrows, ncols;
static DCELL *in_row;
static CELL *old_x_row, *old_y_row;
@@ -52,6 +53,16 @@
return abs(dx) + abs(dy);
}
+static double geodesic_distance(int x1, int y1, int x2, int y2)
+{
+ double lat1 = G_row_to_northing(y1 + 0.5, &window);
+ double lat2 = G_row_to_northing(y2 + 0.5, &window);
+ double lon1 = G_col_to_easting(x1 + 0.5, &window);
+ double lon2 = G_col_to_easting(x2 + 0.5, &window);
+
+ return G_geodesic_distance(lon1, lat1, lon2, lat2);
+}
+
void swap_rows(void)
{
CELL *temp;
@@ -69,7 +80,8 @@
old_val_row = new_val_row;
new_val_row = dtemp;
}
-static void check(int col, int dx, int dy)
+
+static void check(int row, int col, int dx, int dy)
{
const CELL *xrow = dy ? old_x_row : new_x_row;
const CELL *yrow = dy ? old_y_row : new_y_row;
@@ -92,7 +104,9 @@
x = xrow[col + dx] + dx;
y = yrow[col + dx] + dy;
v = vrow[col + dx];
- d = (*distance) (xres * x, yres * y);
+ d = distance
+ ? (*distance) (xres * x, yres * y)
+ : geodesic_distance(col, row, col + x, row + y);
if (!G_is_d_null_value(&dist_row[col]) && dist_row[col] < d)
return;
@@ -110,6 +124,10 @@
{
struct Option *in, *dist, *val, *met;
} opt;
+ struct
+ {
+ struct Flag *m;
+ } flag;
const char *in_name;
const char *dist_name;
const char *val_name;
@@ -122,7 +140,7 @@
struct FPRange range;
DCELL min, max;
DCELL *out_row;
- struct Cell_head window;
+ double scale = 1.0;
G_gisinit(argv[0]);
@@ -136,23 +154,25 @@
opt.dist = G_define_standard_option(G_OPT_R_OUTPUT);
opt.dist->key = "distance";
opt.dist->required = NO;
- opt.dist->description = _("Name for distance output raster map");
- opt.dist->guisection = _("Output");
+ opt.dist->description = _("Name for distance output map");
opt.val = G_define_standard_option(G_OPT_R_OUTPUT);
opt.val->key = "value";
opt.val->required = NO;
- opt.val->description = _("Name for value output raster map");
- opt.val->guisection = _("Output");
+ opt.val->description = _("Name for value output map");
opt.met = G_define_option();
opt.met->key = "metric";
opt.met->type = TYPE_STRING;
opt.met->required = NO;
opt.met->description = _("Metric");
- opt.met->options = "euclidean,squared,maximum,manhattan";
+ opt.met->options = "euclidean,squared,maximum,manhattan,geodesic";
opt.met->answer = "euclidean";
+ flag.m = G_define_flag();
+ flag.m->key = 'm';
+ flag.m->description = _("Output distances in meters instead of map units");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -163,6 +183,8 @@
if (!dist_name && !val_name)
G_fatal_error(_("At least one of distance= and value= must be given"));
+ G_get_window(&window);
+
if (strcmp(opt.met->answer, "euclidean") == 0)
distance = &distance_euclidean_squared;
else if (strcmp(opt.met->answer, "squared") == 0)
@@ -171,6 +193,14 @@
distance = &distance_maximum;
else if (strcmp(opt.met->answer, "manhattan") == 0)
distance = &distance_manhattan;
+ else if (strcmp(opt.met->answer, "geodesic") == 0) {
+ double a, e2;
+ if (window.proj != PROJECTION_LL)
+ G_fatal_error(_("metric=geodesic is only valid for lat/lon"));
+ distance = NULL;
+ G_get_ellipsoid_parameters(&a, &e2);
+ G_begin_geodesic_distance(a, e2);
+ }
else
G_fatal_error(_("Unknown metric: [%s]."), opt.met->answer);
@@ -178,6 +208,12 @@
if (in_fd < 0)
G_fatal_error(_("Unable to open raster map <%s>"), in_name);
+ if (flag.m->answer) {
+ scale = G_database_units_to_meters_factor();
+ if (strcmp(opt.met->answer, "squared") == 0)
+ scale *= scale;
+ }
+
if (dist_name) {
dist_fd = G_open_raster_new(dist_name, DCELL_TYPE);
if (dist_fd < 0)
@@ -195,8 +231,6 @@
if (temp_fd < 0)
G_fatal_error(_("Unable to create temporary file <%s>"), temp_name);
- G_get_window(&window);
-
nrows = window.rows;
ncols = window.cols;
xres = window.ew_res;
@@ -222,7 +256,7 @@
G_set_c_null_value(old_x_row, ncols);
G_set_c_null_value(old_y_row, ncols);
- G_important_message(_("Reading input data..."));
+ G_message(_("Reading raster map <%s>..."), opt.in->answer);
for (row = 0; row < nrows; row++) {
int irow = nrows - 1 - row;
@@ -244,15 +278,15 @@
}
for (col = 0; col < ncols; col++)
- check(col, -1, 0);
+ check(irow, col, -1, 0);
for (col = ncols - 1; col >= 0; col--)
- check(col, 1, 0);
+ check(irow, col, 1, 0);
for (col = 0; col < ncols; col++) {
- check(col, -1, 1);
- check(col, 0, 1);
- check(col, 1, 1);
+ check(irow, col, -1, 1);
+ check(irow, col, 0, 1);
+ check(irow, col, 1, 1);
}
write(temp_fd, new_x_row, ncols * sizeof(CELL));
@@ -262,14 +296,15 @@
swap_rows();
}
- G_percent(1, 1, 1);
+ G_percent(row, nrows, 2);
+
G_close_cell(in_fd);
G_set_c_null_value(old_x_row, ncols);
G_set_c_null_value(old_y_row, ncols);
- G_important_message(_("Writing output data..."));
+ G_message(_("Writing output raster maps..."), opt.in->answer);
for (row = 0; row < nrows; row++) {
int irow = nrows - 1 - row;
off_t offset =
@@ -285,23 +320,27 @@
read(temp_fd, new_val_row, ncols * sizeof(DCELL));
for (col = 0; col < ncols; col++) {
- check(col, -1, -1);
- check(col, 0, -1);
- check(col, 1, -1);
+ check(row, col, -1, -1);
+ check(row, col, 0, -1);
+ check(row, col, 1, -1);
}
for (col = 0; col < ncols; col++)
- check(col, -1, 0);
+ check(row, col, -1, 0);
for (col = ncols - 1; col >= 0; col--)
- check(col, 1, 0);
+ check(row, col, 1, 0);
if (dist_name) {
if (out_row != dist_row)
for (col = 0; col < ncols; col++)
out_row[col] = sqrt(dist_row[col]);
- G_put_d_raster_row(dist_fd, out_row);
+ if (scale != 1.0)
+ for (col = 0; col < ncols; col++)
+ out_row[col] *= scale;
+
+ G_put_d_raster_row(dist_fd, out_row);
}
if (val_name)
@@ -310,7 +349,7 @@
swap_rows();
}
- G_percent(1, 1, 1);
+ G_percent(row, nrows, 2);
close(temp_fd);
remove(temp_name);
More information about the grass-commit
mailing list