[GRASS-SVN] r34429 - in grass/trunk: raster/r.grow.distance
scripts/r.buffer
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 21 10:51:09 EST 2008
Author: glynn
Date: 2008-11-21 10:51:09 -0500 (Fri, 21 Nov 2008)
New Revision: 34429
Modified:
grass/trunk/raster/r.grow.distance/main.c
grass/trunk/scripts/r.buffer/r.buffer.py
Log:
Add metric=geodesic to r.grow.distance
Use metric=geodesic for r.buffer in lat/lon
Modified: grass/trunk/raster/r.grow.distance/main.c
===================================================================
--- grass/trunk/raster/r.grow.distance/main.c 2008-11-21 15:09:33 UTC (rev 34428)
+++ grass/trunk/raster/r.grow.distance/main.c 2008-11-21 15:51:09 UTC (rev 34429)
@@ -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]);
@@ -148,9 +166,13 @@
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);
@@ -161,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)
@@ -169,9 +193,20 @@
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);
+ if (flag.m->answer)
+ scale = G_database_units_to_meters_factor();
+
in_fd = G_open_cell_old(in_name, "");
if (in_fd < 0)
G_fatal_error(_("Unable to open raster map <%s>"), in_name);
@@ -193,8 +228,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;
@@ -241,15 +274,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));
@@ -282,16 +315,26 @@
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(row, col, -1, 0);
+
+ for (col = ncols - 1; col >= 0; col--)
+ 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]);
+ if (scale != 1.0)
+ for (col = 0; col < ncols; col++)
+ out_row[col] *= scale;
+
G_put_d_raster_row(dist_fd, out_row);
}
Modified: grass/trunk/scripts/r.buffer/r.buffer.py
===================================================================
--- grass/trunk/scripts/r.buffer/r.buffer.py 2008-11-21 15:09:33 UTC (rev 34428)
+++ grass/trunk/scripts/r.buffer/r.buffer.py 2008-11-21 15:51:09 UTC (rev 34429)
@@ -99,7 +99,14 @@
distances1 = [scale * float(d) for d in distances]
distances2 = [d * d for d in distances1]
- grass.run_command('r.grow.distance', input = input, metric = 'squared',
+ s = grass.read_command("g.proj", flags='j')
+ kv = grass.parse_key_val(s)
+ if kv['+proj'] == 'longlat':
+ metric = 'geodesic'
+ else:
+ metric = 'squared'
+
+ grass.run_command('r.grow.distance', input = input, metric = metric,
distance = temp_dist)
if zero:
More information about the grass-commit
mailing list