r.profile: feet support (merge from trunk, https://trac.osgeo.org/grass/changeset/62035,r62046,r62060)

git-svn-id: https://svn.osgeo.org/grass/grass/branches/releasebranch_7_0@62177 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))
 

include/gis.h

@@ -74,6 +74,7 @@ static const char *GRASS_copyright __attribute__ ((unused))
 #define U_FEET		6
 #define U_RADIANS	7
 #define U_DEGREES	8
+#define U_USFEET	9
 /* Temporal units from the datetime library */
 #define U_YEARS         DATETIME_YEAR   
 #define U_MONTHS        DATETIME_MONTH  

lib/gis/proj2.c

@@ -23,6 +23,7 @@
    - U_UNKNOWN (XY)
    - U_METERS  (UTM)
    - U_FEET    (SP)
+   - U_USFEET (a few SP)
    - U_DEGREES (LL)
    
   \return units code (see gis.h)
@@ -36,7 +37,7 @@ int G__projection_units(int n)
     case PROJECTION_UTM:
 	return U_METERS;
     case PROJECTION_SP:
-	return U_FEET;
+	return U_FEET; /* TODO: what if U_USFEET as in CA and NC ? */
     case PROJECTION_LL:
 	return U_DEGREES;
     default:

lib/gis/proj3.c

@@ -73,6 +73,8 @@ const char *G_database_unit_name(int plural)
 	    units = U_MILES;
 	else if (strcasecmp(name, "foot") == 0 || strcasecmp(name, "feet") == 0)
 	    units = U_FEET;
+	else if (strcasecmp(name, "foot_us") == 0 || strcasecmp(name, "foot_uss") == 0)
+	    units = U_USFEET;
 	else if (strcasecmp(name, "degree") == 0 || strcasecmp(name, "degrees") == 0)
 	    units = U_DEGREES;
 	else
@@ -127,6 +129,7 @@ double G_database_units_to_meters_factor(void)
     double factor;
     int n;
 
+    /* TODO: sync with definitions in ../proj/units.table */
     static const struct
     {
 	char *unit;
@@ -134,7 +137,8 @@ double G_database_units_to_meters_factor(void)
     } table[] = {
 	{"unit",  1.0},
 	{"meter", 1.0},
-	{"foot", .3048},
+	{"foot",  .3048},
+	{"foot_us", 1200/3937.},
 	{"inch", .0254},
 	{NULL, 0.0}
     };

lib/gis/units.c

@@ -26,6 +26,7 @@
    - U_KILOMETERS
    - U_MILES
    - U_FEET
+   - U_USFEET
 
   Returns a factor which converts meters to units (by multiplication).
  
@@ -52,6 +53,10 @@ double G_meters_to_units_factor(int units)
 	return 3.28083989501312;	        /*  1 / (0.0254 * 12)    */
 	break;
 	
+    case U_USFEET:
+	return 3.28083333333333;       	        /*  1 / (1200/3937)    */
+	break;
+	
     default:
 	return 1.0;
 	break;
@@ -70,6 +75,7 @@ double G_meters_to_units_factor(int units)
    - U_HECTARES
    - U_MILES
    - U_FEET
+   - U_USFEET
 
   Returns a factor which converts square meters to square units (by
   multiplication).
@@ -102,7 +108,11 @@ double G_meters_to_units_factor_sq(int units)
 	break;
 	
     case U_FEET:
-	return 10.7639104167097;	        /*  1 / (0.0254 * 12)^2  */
+	return 10.7639104167097;	/*  1 / (0.0254 * 12)^2  */
+	break;
+	
+    case U_USFEET:
+	return 10.7638673611111;       	/*  1 / (1200/3937)^2    */
 	break;
 	
     default:
@@ -135,6 +145,8 @@ int G_is_units_type_spatial(int units)
         return 1;
     case U_FEET:
         return 1;
+    case U_USFEET:
+        return 1;
     case U_RADIANS:
         return 1;
     case U_DEGREES:
@@ -179,6 +191,7 @@ int G_is_units_type_temporal(int units)
    - U_HECTARES
    - U_MILES
    - U_FEET
+   - U_USFEET
 
   \param units units code
   \param plural plural form if true
@@ -241,6 +254,13 @@ const char *G_get_units_name(int units, int plural, int square)
 	    return plural ? _("feet") : _("foot");
 	break;
 
+    case U_USFEET:
+	if (square)
+	    return plural ? _("square US feet") : _("square US foot");
+	else
+	    return plural ? _("US feet") : _("US foot");
+	break;
+
     case U_DEGREES:
 	if (square)
 	    return plural ? _("square degrees") : _("square degree");
@@ -286,6 +306,7 @@ const char *G_get_units_name(int units, int plural, int square)
    - U_HECTARES
    - U_MILES
    - U_FEET
+   - U_USFEET
    - ...
    - U_YEARS
    - ...
@@ -319,11 +340,14 @@ int G_units(const char *units_name)
     else if (strcasecmp(units_name, "foot") == 0 ||
 	     strcasecmp(units_name, "feet") == 0)
 	return U_FEET;
+    else if (strcasecmp(units_name, "foot_us") == 0 ||
+	     strcasecmp(units_name, "foot_uss") == 0)
+	return U_USFEET;
     else if (strcasecmp(units_name, "degree") == 0 ||
 	     strcasecmp(units_name, "degrees") == 0)
 	return U_DEGREES;
     else if (strcasecmp(units_name, "year") == 0 ||
-	strcasecmp(units_name, "years") == 0)
+	     strcasecmp(units_name, "years") == 0)
 	return U_YEARS;
     else if (strcasecmp(units_name, "month") == 0 ||
 	     strcasecmp(units_name, "months") == 0)

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,9 @@ static double dist, e, n;
 int main(int argc, char *argv[])
 {
     char *name, *outfile;
+    const char *unit;
+    int unit_id;
+    double factor;
     int fd, projection;
     FILE *fp, *coor_fp;
     double res;
@@ -37,8 +40,8 @@ int main(int argc, char *argv[])
     struct Cell_head window;
     struct
     {
-	struct Option *opt1, *profile, *res, *output, *null_str, *coord_file;
-	struct Flag *g, *c;
+	struct Option *opt1, *profile, *res, *output, *null_str, *coord_file, *units;
+	struct Flag *g, *c, *m;
     }
     parm;
     struct GModule *module;
@@ -98,11 +101,15 @@ int main(int argc, char *argv[])
     parm.c->description =
 	_("Output RRR:GGG:BBB color values for each profile point");
 
+    parm.units = G_define_standard_option(G_OPT_M_UNITS);
+    parm.units->options = "meters,kilometers,feet,miles";
+    parm.units->label = parm.units->description;
+    parm.units->description = _("If units are not specified, current location units are used. "
+                                "Meters are used by default in geographic (latlon) locations.");
 
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-
     clr = 0;
     if (parm.c->answer)
 	clr = 1;		/* color output */
@@ -116,18 +123,36 @@ int main(int argc, char *argv[])
 
     G_get_window(&window);
     projection = G_projection();
+
+    /* get conversion factor and units name */
+    if (parm.units->answer) {
+        unit_id = G_units(parm.units->answer);
+        factor = 1. / G_meters_to_units_factor(unit_id);
+        unit = G_get_units_name(unit_id, 1, 0);
+    }
+    /* keep meters in case of latlon */
+    else if (projection == PROJECTION_LL) {
+        factor = 1;
+        unit = "meters";
+    } 
+    else {
+        /* get conversion factor to current units */
+        unit = G_database_unit_name(1);
+        factor = G_database_units_to_meters_factor();
+    }
+
     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 +183,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 +211,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 +236,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 +250,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 +268,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 +279,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 +289,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 +305,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 +327,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 +346,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 +365,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>Option <b>units</b> enables to set units of the profile length output.
+If the units are not specified, current location units will be used.
+In case of geographic locations (latitude/longitude), meters are used as default unit.
 
 <h2>EXAMPLES</h2>