[GRASS-SVN] r58550 - grass/trunk/raster/r.horizon
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Dec 29 15:09:07 PST 2013
Author: neteler
Date: 2013-12-29 15:09:07 -0800 (Sun, 29 Dec 2013)
New Revision: 58550
Modified:
grass/trunk/raster/r.horizon/TODO
grass/trunk/raster/r.horizon/main.c
grass/trunk/raster/r.horizon/r.horizon.html
Log:
r.horizon: added sanity check for dist parameter; fixed wrong user msg; format single point mode output as CSV; expand parameter names properly; update manual examples to NC
Modified: grass/trunk/raster/r.horizon/TODO
===================================================================
--- grass/trunk/raster/r.horizon/TODO 2013-12-29 22:30:51 UTC (rev 58549)
+++ grass/trunk/raster/r.horizon/TODO 2013-12-29 23:09:07 UTC (rev 58550)
@@ -1,8 +1,5 @@
TODO
-11/2008: Integrate latest GRASS 7 API changes
-
--------
Probably the sun position calculation should be replaced
with
Modified: grass/trunk/raster/r.horizon/main.c
===================================================================
--- grass/trunk/raster/r.horizon/main.c 2013-12-29 22:30:51 UTC (rev 58549)
+++ grass/trunk/raster/r.horizon/main.c 2013-12-29 23:09:07 UTC (rev 58550)
@@ -24,8 +24,9 @@
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ * along with this program; if not, write to the
+ * Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <stdio.h>
@@ -57,16 +58,14 @@
FILE *fw;
-const double rad2deg = 180. / M_PI;
-const double deg2rad = M_PI / 180.;
const double pihalf = M_PI * 0.5;
-const double twopi = 2. * M_PI;
+const double twopi = M_PI * 2.;
const double invEarth = 1. / EARTHRADIUS;
-
+const double deg2rad = M_PI / 180.;
+const double rad2deg = 180. / M_PI;
const double minAngle = DEG;
const char *elevin;
-const char *latin = NULL;
const char *horizon = NULL;
const char *mapset = NULL;
const char *per;
@@ -141,7 +140,8 @@
int ll_correction = FALSE;
double coslatsq;
-/* use G_distance() instead ??!?! */
+/* why not use G_distance() here which switches to geodesic/great
+ circle distance as needed? */
double distance(double x1, double x2, double y1, double y2)
{
if (ll_correction) {
@@ -189,7 +189,7 @@
" counterclockwise with east=0, north=90 etc. The output is the horizon height in radians.");
parm.elevin = G_define_option();
- parm.elevin->key = "elevin";
+ parm.elevin->key = "elev_in";
parm.elevin->type = TYPE_STRING;
parm.elevin->required = YES;
parm.elevin->gisprompt = "old,cell,raster";
@@ -207,7 +207,7 @@
parm.direction->guisection = _("Input options");
parm.step = G_define_option();
- parm.step->key = "horizonstep";
+ parm.step->key = "horizon_step";
parm.step->type = TYPE_DOUBLE;
parm.step->required = NO;
parm.step->description =
@@ -273,7 +273,7 @@
parm.coord = G_define_option();
- parm.coord->key = "coord";
+ parm.coord->key = "coordinate";
parm.coord->type = TYPE_DOUBLE;
parm.coord->key_desc = "east,north";
parm.coord->required = NO;
@@ -282,7 +282,7 @@
parm.coord->guisection = _("Output options");
parm.dist = G_define_option();
- parm.dist->key = "dist";
+ parm.dist->key = "distance";
parm.dist->type = TYPE_DOUBLE;
parm.dist->answer = DIST;
parm.dist->required = NO;
@@ -341,7 +341,7 @@
setMode(SINGLE_POINT);
if (sscanf(parm.coord->answer, "%lf,%lf", &xcoord, &ycoord) != 2) {
G_fatal_error(
- _("Can't read the coordinates from the \"coord\" option."));
+ _("Can't read the coordinates from the \"coordinate\" option."));
}
/* Transform the coordinates to row/column */
@@ -423,6 +423,7 @@
sscanf(parm.dist->answer, "%lf", &dist);
+ if (dist < 0.5 || dist > 1.5 ) G_fatal_error(_("The distance value must be 0.5-1.5. Aborting."));
stepxy = dist * 0.5 * (stepx + stepy);
distxy = dist;
@@ -469,8 +470,7 @@
if ((in_proj_info = G_get_projinfo()) == NULL)
G_fatal_error(
- _("Can't get projection info of current location: "
- "please set latitude via 'lat' or 'latin' option!"));
+ _("Can't get projection info of current location"));
if ((in_unit_info = G_get_projunits()) == NULL)
G_fatal_error(_("Can't get projection units of current location"));
@@ -495,6 +495,7 @@
INPUT();
+ G_debug(3, "calculate() starts...");
calculate(xcoord, ycoord, (int)(ebufferZone / stepx),
(int)(wbufferZone / stepx), (int)(sbufferZone / stepy),
(int)(nbufferZone / stepy));
@@ -826,7 +827,7 @@
else if (printangle >= 360.)
printangle -= 360;
- G_message("%lf, %lf", printangle, shadow_angle);
+ G_message("%lf,%lf", printangle, shadow_angle);
angle += dfr_rad;
@@ -1212,6 +1213,7 @@
}
}
+ G_debug(3, "OUTGR() starts...");
OUTGR(cellhd.rows, cellhd.cols);
/* empty array */
Modified: grass/trunk/raster/r.horizon/r.horizon.html
===================================================================
--- grass/trunk/raster/r.horizon/r.horizon.html 2013-12-29 22:30:51 UTC (rev 58549)
+++ grass/trunk/raster/r.horizon/r.horizon.html 2013-12-29 23:09:07 UTC (rev 58550)
@@ -9,11 +9,12 @@
heights in the specified directions from the given point. The results are
written to the stdout.
<li> raster: in this case the output is
-one or more rasters, with each point in a raster giving the horizon
+one or more raster maps, with each point in a raster giving the horizon
height in a specific direction. One raster is created for each direction.
</ul>
-<p>The directions are given as azimuthal angles (in degrees), with
+<p>
+The directions are given as azimuthal angles (in degrees), with
the angle starting with 0 towards East and moving counterclockwise
(North is 90, etc.). The calculation takes into account the actual
projection, so the angles are corrected for direction distortions
@@ -34,37 +35,40 @@
<h3>Input parameters:</h3>
-<p>The <I>elevin</I> parameter is an input elevation raster map. If
+<p>The <i>elev_in</i> parameter is an input elevation raster map. If
the buffer options are used (see below), this raster should extend
over the area that accommodate the presently defined region plus
defined buffer zones.
-<p>The <I>horizonstep</I> parameter gives the angle step (in degrees)
+<p>The <i>horizon_step</i> parameter gives the angle step (in degrees)
between successive azimuthal directions for the calculation of the
-horizon. Thus, a value of 5 for the <I>horizonstep</I> will give a total of
-360/5=72 directions (72 rasters if used in the raster mode).
+horizon. Thus, a value of 5 for the <i>horizon_step</i> will give a total of
+360/5=72 directions (72 raster maps if used in the raster mode).
-<p>The <I>direction</I> parameter gives the initial direction of the
+<p>The <i>direction</i> parameter gives the initial direction of the
first output. This parameter acts as an direction angle offset. For
example, if you want to get horizon angles for directions 45 and 225
-degrees, the <I>direction</I> should be set to 45 and <I>horizonstep</I> to
+degrees, the <i>direction</i> should be set to 45 and <i>horizon_step</i> to
180. If you only want one single direction, use this parameter to
-specify desired direction of horizon angle, and set the <I>horizonstep</I>
-size to 0 degrees.
+specify desired direction of horizon angle, and set the <i>horizon_step</i>
+size to 0 degrees. Otherwise all angles for a given starting <i>direction</i>
+with step of <i>horizon_step</i> are calculated.
-<p>The <I>dist </I>controls the sampling distance step size for the
+<p>The <i>distance</i> controls the sampling distance step size for the
search for horizon along the line of sight. The default value is 1.0
meaning that the step size will be taken from the raster resolution.
Setting the value below 1.0 might slightly improve results for
directions apart from the cardinal ones, but increasing the
processing load of the search algorithm.
-<p>The <I>maxdistance</I> value gives a maximum distance to move away
+<p>The <i>maxdistance</i> value gives a maximum distance to move away
from the origin along the line of sight in order to search for the
horizon height. The smaller this value the faster the calculation but
the higher the risk that you may miss a terrain feature that can
-contribute significantly to the horizon outline.
-<p>The <I>coord</I> parameter takes a pair of easting-northing values
+contribute significantly to the horizon outline. Note that a viewshed
+can be calculated with <em>r.viewshed</em>.
+
+<p>The <i>coordinate</i> parameter takes a pair of easting-northing values
in the current coordinate system and calculates the values of angular
height of the horizon around this point. To achieve the
consistency of the results, the point coordinate is
@@ -76,16 +80,16 @@
have contributed to the horizon, because these features are outside
the region. There are to options how to set the size of the buffer
that is used to increase the area of the horizon analysis. The
-<I>bufferzone</I> parameter allows you to specify the same size of
-buffer for all cardinal directions and the parameters <I>e_buff</I>,
-<I>n_buff</I>, <I>s_buff</I>, and <I>w_buff</I> allow you to specify
+<i>bufferzone</i> parameter allows you to specify the same size of
+buffer for all cardinal directions and the parameters <i>e_buff</i>,
+<i>n_buff</i>, <i>s_buff</i>, and <i>w_buff</i> allow you to specify
a buffer size individually for each of the four directions. The
buffer parameters influence only size of the read elevation map,
while the analysis in the raster mode will be done only for the
area specified by the current region definition.
-<p>The <I>horizon </I>parameter gives the prefix of the output
+<p>The <i>horizon </i>parameter gives the prefix of the output
horizon raster maps. The raster name of each horizon direction
-raster will be constructed as <I>horizon_</I>NNN , where NNN counts
+raster will be constructed as <i>horizon_</i>NNN , where NNN counts
upwards from 0 to total number of directions. If you use <b>r.horizon</b>
in the single point mode this option will be ignored.
@@ -106,30 +110,37 @@
allow the line of sight to pass just above the terrain at that point.
This is continued until the line of sight reaches a height that is
higher than any point in the region or until it reaches the border of
-the region (see also the <I>bufferzone,e_buff</I>, <I>n_buff</I>,
-<I>s_buff</I>, and <I>w_buff</I>). The the number of lines of sight (azimuth
-directions) is determined from the <I>direction</I> and
-<I>horizonstep</I> parameters. The method takes into account the curvature
+the region (see also the <i>bufferzone,e_buff</i>, <i>n_buff</i>,
+<i>s_buff</i>, and <i>w_buff</i>). The the number of lines of sight (azimuth
+directions) is determined from the <i>direction</i> and
+<i>horizon_step</i> parameters. The method takes into account the curvature
of the Earth whereby remote features will seem to be lower than they
actually are. It also accounts for the changes of angles towards
cardinal directions caused by the projection (see above).
-<h2>EXAMPLE</h2>
+<h2>EXAMPLES</h2>
-Single point mode:
+The examples are intended for the North Carolina sample dataset:
+<p>
+Single point mode (output of horizon angles CCW from East):
<div class="code"><pre>
-r.horizon elevin=DEM horizonstep=30 direction=15 bufferzone=200 \
- coord=47.302,7.365 dist=0.7 > horizon.out
+# determine horizon angle in 225 degree direction:
+r.horizon elev_in=elevation direction=215 horizon_step=0 bufferzone=200 \
+ coordinate=638871.6,223384.4 maxdistance=5000
+
+# determine horizon values starting at 90 deg (North), step size of 5 deg:
+r.horizon elev_in=elevation direction=90 horizon_step=5 bufferzone=200 \
+ coordinate=638871.6,223384.4 maxdistance=5000
</pre></div>
-Raster map mode (for r.sun):
+Raster map mode (output maps "horangle*" become input for <em>r.sun</em>):
<div class="code"><pre>
# we put a bufferzone of 10% of maxdistance around the study area
-r.horizon elevin=DEM horizonstep=30 bufferzone=200 horizon=horangle \
- dist=0.7 maxdistance=2000
+r.horizon elev_in=elevation horizon_step=30 bufferzone=200 horizon=horangle \
+ maxdistance=5000
</pre></div>
@@ -138,7 +149,8 @@
<em>
<a href="r.sun.html">r.sun</a>,
<a href="r.sunmask.html">r.sunmask</a>,
-<a href="r.los.html">r.los</a></em>
+<a href="r.los.html">r.los</a>,
+<a href="r.viewshed.html">r.viewshed</a></em>
<h2>REFERENCES</h2>
@@ -181,4 +193,4 @@
<a href="mailto:Marcel.Suri at jrc.it">Marcel.Suri at jrc.it</a>
</ADDRESS>
-<p><I>Last changed: $Date$</I>
+<p><i>Last changed: $Date$</i>
More information about the grass-commit
mailing list