Add metric=geodesic to r.grow.distance

Use metric=geodesic for r.buffer in lat/lon


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

raster/r.grow.distance/main.c

@@ -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 @@ static double distance_manhattan(double dx, double dy)
     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 @@ void swap_rows(void)
     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 @@ static void check(int col, int dx, int dy)
     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 @@ int main(int argc, char **argv)
     {
 	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 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
 	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 @@ int main(int argc, char **argv)
     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 @@ int main(int argc, char **argv)
 	    }
 
 	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 @@ int main(int argc, char **argv)
 	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);
 	}
 

scripts/r.buffer/r.buffer.py

@@ -99,7 +99,14 @@ def main():
     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: