r.profile: use location units, revert gui profile tool changes, add r.profile test, see https://trac.osgeo.org/grass/ticket/2417

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

gui/wxpython/wxplot/profile.py

@@ -44,21 +44,6 @@ from wxplot.dialogs    import ProfileRasterDialog, PlotStatsFrame
 from core.gcmd         import RunCommand, GWarning, GError, GMessage
 
 
-def get_unit_conversion():
-    """Get conversion factor to convert output from r.profile to current units.
-    Module r.profile always outputs distance in meters.
-
-    .. todo:: create a wrapper for g.proj in grass.core    
-    """
-    proj = RunCommand(prog='g.proj', read=True, flags='g').strip()
-    if not proj:
-        # TODO put some warning here
-        return 1.
-
-    proj = grass.parse_key_val(proj, sep='=')
-    return 1 / float(proj.get('meters', 1))
-
-
 class ProfileFrame(BasePlotFrame):
     """Mainframe for displaying profile of one or more raster maps. Uses wx.lib.plot.
     """
@@ -74,8 +59,8 @@ class ProfileFrame(BasePlotFrame):
         if sys.platform != 'darwin':
             self.SetToolBar(self.toolbar)
         self.SetTitle(_("GRASS Profile Analysis Tool"))
-        self._units = units
-        self._conversion = get_unit_conversion()
+        # in case of degrees, we want to use meters
+        self._units = units if 'degree' not in units else 'meters'
 
         #
         # Init variables
@@ -275,8 +260,7 @@ class ProfileFrame(BasePlotFrame):
             if dist == None or dist == '' or dist == 'nan' or \
                 elev == None or elev == '' or elev == 'nan':
                 continue
-            # remove conversion if one day r.profile starts to output distance in projection units!
-            dist = self._conversion * float(dist)
+            dist = float(dist)
             elev = float(elev)
             datalist.append((dist,elev))
 

raster/r.profile/local_proto.h

@@ -6,8 +6,8 @@
 #include <grass/raster.h>
 
 /* main.c */
-int do_profile(double, double, double, double, char *, int, double, int, int,
-	       FILE *, char *);
+int do_profile(double, double, double, double, int, double, int, int,
+	       FILE *, char *, const char *, double);
 
 /* read_rast.c */
 int read_rast(double, double, double, int, int, RASTER_MAP_TYPE, FILE *,

raster/r.profile/main.c

@@ -23,6 +23,8 @@ static double dist, e, n;
 int main(int argc, char *argv[])
 {
     char *name, *outfile;
+    const char *unit;
+    double factor;
     int fd, projection;
     FILE *fp, *coor_fp;
     double res;
@@ -38,7 +40,7 @@ int main(int argc, char *argv[])
     struct
     {
 	struct Option *opt1, *profile, *res, *output, *null_str, *coord_file;
-	struct Flag *g, *c;
+	struct Flag *g, *c, *m;
     }
     parm;
     struct GModule *module;
@@ -98,6 +100,11 @@ int main(int argc, char *argv[])
     parm.c->description =
 	_("Output RRR:GGG:BBB color values for each profile point");
 
+    parm.m = G_define_flag();
+    parm.m->key = 'm';
+    parm.m->label = _("Use meters instead of current location units");
+    parm.m->description =
+            _("Meters are used always in latlon locations");
 
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
@@ -116,18 +123,29 @@ int main(int argc, char *argv[])
 
     G_get_window(&window);
     projection = G_projection();
+
+    /* get conversion factor to current units */
+    unit = G_database_unit_name(1);
+    factor = G_database_units_to_meters_factor();
+    /* keep meters in case of latlon */
+    if (parm.m->answer || projection == 3)
+    {
+        factor = 1;
+        unit = "meters";
+    }
+
     if (parm.res->answer) {
 	res = atof(parm.res->answer);
 	/* Catch bad resolution ? */
 	if (res <= 0)
-	    G_fatal_error(_("Illegal resolution! [%g]"), res);
+	    G_fatal_error(_("Illegal resolution %g [%s]"), res / factor, unit);
     }
     else {
 	/* Do average of EW and NS res */
 	res = (window.ew_res + window.ns_res) / 2;
     }
 
-    G_message(_("Using resolution [%g]"), res);
+    G_message(_("Using resolution: %g [%s]"), res / factor, unit);
 
     G_begin_distance_calculations();
 
@@ -158,14 +176,14 @@ int main(int argc, char *argv[])
     /* Done with file */
 
     /* Show message giving output format */
-    G_message(_("Output Format:"));
+    G_message(_("Output columns:"));
     if (coords == 1)
 	sprintf(formatbuff,
-		_("[Easting] [Northing] [Along Track Dist.(m)] [Elevation]"));
+		_("Easting, Northing, Along track dist. [%s], Elevation"), unit);
     else
-	sprintf(formatbuff, _("[Along Track Dist.(m)] [Elevation]"));
+	sprintf(formatbuff, _("Along track dist. [%s], Elevation"), unit);
     if (clr)
-	strcat(formatbuff, _(" [RGB Color]"));
+	strcat(formatbuff, _(" RGB color"));
     G_message(formatbuff);
 
     /* Get Profile Start Coords */
@@ -186,8 +204,8 @@ int main(int argc, char *argv[])
 		G_fatal_error(_("Invalid coordinates %s %s"), ebuf, nbuf);
 
 	    if (havefirst)
-		do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type,
-			   fp, null_string);
+		do_profile(e1, e2, n1, n2, coords, res, fd, data_type,
+			   fp, null_string, unit, factor);
 	    e1 = e2;
 	    n1 = n2;
 	    havefirst = TRUE;
@@ -211,8 +229,8 @@ int main(int argc, char *argv[])
 	    n2 = n1;
 
 	    /* Get profile info */
-	    do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type, fp,
-		       null_string);
+	    do_profile(e1, e2, n1, n2, coords, res, fd, data_type, fp,
+		       null_string, unit, factor);
 	}
 	else {
 	    for (i = 0; i <= k - 2; i += 2) {
@@ -225,8 +243,8 @@ int main(int argc, char *argv[])
 				G_projection());
 
 		/* Get profile info */
-		do_profile(e1, e2, n1, n2, name, coords, res, fd, data_type,
-			   fp, null_string);
+		do_profile(e1, e2, n1, n2, coords, res, fd, data_type,
+			   fp, null_string, unit, factor);
 
 	    }
 	}
@@ -243,9 +261,9 @@ int main(int argc, char *argv[])
 
 /* Calculate the Profile Now */
 /* Establish parameters */
-int do_profile(double e1, double e2, double n1, double n2, char *name,
+int do_profile(double e1, double e2, double n1, double n2,
 	       int coords, double res, int fd, int data_type, FILE * fp,
-	       char *null_string)
+	       char *null_string, const char *unit, double factor)
 {
     float rows, cols, LEN;
     double Y, X, AZI;
@@ -254,7 +272,7 @@ int do_profile(double e1, double e2, double n1, double n2, char *name,
     rows = n1 - n2;
 
     LEN = G_distance(e1, n1, e2, n2);
-    G_message(_("Approx. transect length [%f] m"), LEN);
+    G_message(_("Approx. transect length: %f [%s]"), LEN / factor, unit);
 
     if (!G_point_in_region(e2, n2))
 	G_warning(_("Endpoint coordinates are outside of current region settings"));
@@ -264,7 +282,7 @@ int do_profile(double e1, double e2, double n1, double n2, char *name,
 	/* Special case for no movement */
 	e = e1;
 	n = n1;
-	read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+	read_rast(e, n, dist / factor, fd, coords, data_type, fp, null_string);
     }
 
     if (rows >= 0 && cols < 0) {
@@ -280,7 +298,7 @@ int do_profile(double e1, double e2, double n1, double n2, char *name,
 	    dist -= G_distance(e, n, e1, n1);
 	}
 	for (e = e1, n = n1; e < e2 || n > n2; e += X, n -= Y) {
-	    read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+	    read_rast(e, n, dist / factor, fd, coords, data_type, fp, null_string);
 	    /* d+=res; */
 	    dist += G_distance(e - X, n + Y, e, n);
 	}
@@ -302,7 +320,7 @@ int do_profile(double e1, double e2, double n1, double n2, char *name,
 	     */
 	}
 	for (e = e1, n = n1; e < e2 || n < n2; e += X, n += Y) {
-	    read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+	    read_rast(e, n, dist / factor, fd, coords, data_type, fp, null_string);
 	    /* d+=res; */
 	    dist += G_distance(e - X, n - Y, e, n);
 	}
@@ -321,7 +339,7 @@ int do_profile(double e1, double e2, double n1, double n2, char *name,
 	    dist -= G_distance(e, n, e1, n1);
 	}
 	for (e = e1, n = n1; e > e2 || n > n2; e -= X, n -= Y) {
-	    read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+	    read_rast(e, n, dist / factor, fd, coords, data_type, fp, null_string);
 	    /* d+=res; */
 	    dist += G_distance(e + X, n + Y, e, n);
 	}
@@ -340,7 +358,7 @@ int do_profile(double e1, double e2, double n1, double n2, char *name,
 	    dist -= G_distance(e, n, e1, n1);
 	}
 	for (e = e1, n = n1; e > e2 || n < n2; e -= X, n += Y) {
-	    read_rast(e, n, dist, fd, coords, data_type, fp, null_string);
+	    read_rast(e, n, dist / factor, fd, coords, data_type, fp, null_string);
 	    /* d+=res; */
 	    dist += G_distance(e + X, n - Y, e, n);
 	}

raster/r.profile/r.profile.html

@@ -2,8 +2,8 @@
 
 This program outputs two or four column (with <b>-g</b>) data to stdout or 
 an ASCII file. The default two column output consists of cumulative profile 
-length (in meters) and raster value. The optional four column output consists 
-of easting, northing, cumlative profile length (m), and raster value. Profile
+length and raster value. The optional four column output consists 
+of easting, northing, cumulative profile length, and raster value. Profile
 end or "turning" points can be set manually with the <b>coordinates</b>
 argument. The profile resolution, or distance between profile
 points, is obtained from the current region resolution, or can be manually
@@ -45,6 +45,9 @@ The 2 column output is compatible with most plotting programs.
 The optional RGB output provides the associated GRASS colour value for
 each profile point.
 
+<p>Profile distance units correspond to the location units.
+Flag <b>-m</b> can be used to force output values in meters.
+In case of latlon locations and locations using meters, this flag does nothing.
 
 <h2>EXAMPLES</h2>
 

raster/r.profile/testsuite/test_profile_ncspm.py

@@ -0,0 +1,127 @@
+from grass.gunittest import TestCase, test
+from grass.gunittest.gmodules import SimpleModule
+import grass.script.core as gcore
+
+# not used yet
+LOCATION = 'nc_spm'
+
+output1 = """
+ 0.000000 88.370453
+ 10.000000 88.397057
+ 20.000000 89.526253
+ 30.000000 89.677551
+ 40.000000 91.297195
+ 50.000000 91.297195
+ 60.000000 92.330658
+ 70.000000 93.069199
+ 80.000000 94.768280
+ 90.000000 95.524551
+ 100.000000 96.770805
+ 110.000000 96.770805
+ 120.000000 97.418869
+"""
+
+output2 = """
+637656.000000 224222.000000 0.000000 88.370453
+637664.540486 224227.201932 10.000000 88.397057
+637673.080972 224232.403865 20.000000 89.526253
+637681.621458 224237.605797 30.000000 89.677551
+637690.161944 224242.807729 40.000000 91.297195
+637698.702430 224248.009662 50.000000 91.297195
+637707.242916 224253.211594 60.000000 92.330658
+637715.783402 224258.413526 70.000000 93.069199
+637724.323887 224263.615459 80.000000 94.768280
+637732.864373 224268.817391 90.000000 95.524551
+637741.404859 224274.019323 100.000000 96.770805
+637749.945345 224279.221256 110.000000 96.770805
+637758.485831 224284.423188 120.000000 97.418869
+"""
+
+output3 = """
+ 0.000000 91.071831
+ 10.000000 91.431198
+ 20.000000 91.746628
+ 30.000000 91.746628
+ 40.000000 91.748047
+ 50.000000 91.872192
+ 60.000000 91.730049
+ 70.000000 91.690292
+ 80.000000 91.341331
+ 86.533231 91.341331
+ 96.533231 91.639000
+ 106.533231 nodata
+ 116.533231 nodata
+ 126.533231 nodata
+ 136.533231 nodata
+ 146.533231 nodata
+ 156.533231 nodata
+ 166.533231 nodata
+ 176.533231 nodata
+ 186.533231 nodata
+ 196.533231 nodata
+ 206.533231 nodata
+ 216.533231 nodata
+"""
+
+output4 = """
+ 0.000000 88.370453
+ 25.000000 89.526253
+ 50.000000 91.297195
+ 75.000000 94.768280
+ 100.000000 96.770805
+ 125.000000 97.646629
+"""
+
+
+class TestProfileNCSPM(TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        gcore.use_temp_region()
+        gcore.run_command('g.region', rast='elevation')
+
+    @classmethod
+    def tearDownClass(cls):
+        gcore.del_temp_region()
+
+    def test_profile_default(self):
+        rprofile = SimpleModule('r.profile', input='elevation',
+                                coordinates=[637656, 224222, 637766, 224289])
+        self.assertModule(rprofile)
+        self.assertMultiLineEqual(rprofile.outputs.stdout.strip(), output1.strip())
+        self.assertIn('128.798294 [meters]', rprofile.outputs.stderr)  # distance
+        self.assertIn('10 [meters]', rprofile.outputs.stderr)  # resolution
+
+    def test_profile_m(self):
+        rprofile = SimpleModule('r.profile', input='elevation', flags='m',
+                                coordinates=[637656, 224222, 637766, 224289])
+        self.assertModule(rprofile)
+        self.assertIn('128.798294 [meters]', rprofile.outputs.stderr)  # distance
+        self.assertIn('10 [meters]', rprofile.outputs.stderr)  # resolution
+
+    def test_profile_resolution(self):
+        rprofile = SimpleModule('r.profile', input='elevation', resolution=25,
+                                coordinates=[637656, 224222, 637766, 224289])
+        self.assertModule(rprofile)
+        self.assertMultiLineEqual(rprofile.outputs.stdout.strip(), output4.strip())
+        self.assertIn('128.798294 [meters]', rprofile.outputs.stderr)  # distance
+        self.assertIn('25 [meters]', rprofile.outputs.stderr)  # resolution
+
+    def test_profile_ne(self):
+        rprofile = SimpleModule('r.profile', input='elevation', flags='g',
+                                coordinates=[637656, 224222, 637766, 224289])
+        self.assertModule(rprofile)
+        self.assertMultiLineEqual(rprofile.outputs.stdout.strip(), output2.strip())
+
+    def test_profile_region(self):
+        rprofile = SimpleModule('r.profile', input='elevation', null='nodata',
+                                coordinates=[644914, 224579, 644986,
+                                             224627, 645091, 224549])
+        self.assertModule(rprofile)
+        self.assertMultiLineEqual(rprofile.outputs.stdout.strip(), output3.strip())
+        self.assertIn("WARNING: Endpoint coordinates are outside of current region settings",
+                      rprofile.outputs.stderr,)
+
+
+if __name__ == '__main__':
+    test()