[GRASS-SVN] r38685 - grass/branches/develbranch_6/raster/r.horizon
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Aug 11 07:58:03 EDT 2009
Author: hamish
Date: 2009-08-11 07:58:02 -0400 (Tue, 11 Aug 2009)
New Revision: 38685
Modified:
grass/branches/develbranch_6/raster/r.horizon/description.html
grass/branches/develbranch_6/raster/r.horizon/main.c
Log:
write out metadata; whitespace/msg cleanup
Modified: grass/branches/develbranch_6/raster/r.horizon/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/description.html 2009-08-11 11:50:41 UTC (rev 38684)
+++ grass/branches/develbranch_6/raster/r.horizon/description.html 2009-08-11 11:58:02 UTC (rev 38685)
@@ -26,12 +26,14 @@
LOCATION with x,y coordinate system, where this correction is
not applied.
+
<H3>Flags:</H3>
<dl>
<dt><B>-d</B>
<dd>Output horizon height in degrees (the default is radians)</dl>
</dl>
+
<H3>Input parameters:</H3>
<P>The <I>elevin</I> parameter is an input elevation raster map. If
the buffer options are used (see below), this raster should extend
@@ -97,6 +99,7 @@
raster coordinates.
</p>
+
<H2>METHOD</H2>
<P>The calculation method is based on the method used in <B>r.sun</B>
to calculate shadows. It starts at a very shallow angle and walks
@@ -114,18 +117,21 @@
cardinal directions caused by the projection (see above).
</P>
+
<H2>EXAMPLE</H2>
Single point mode:
<div class="code"><pre>
-r.horizon elevin=DEM horizonstep=30 direction=15 bufferzone=200 coord=47.302,7.365 dist=0.7 > horizon.out
+r.horizon elevin=DEM horizonstep=30 direction=15 bufferzone=200 \
+ coord=47.302,7.365 dist=0.7 > horizon.out
</pre></div>
Raster map mode (for r.sun):
<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 elevin=DEM horizonstep=30 bufferzone=200 horizon=horangle \
+ dist=0.7 maxdistance=2000
</pre></div>
@@ -133,8 +139,10 @@
<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>
+
<H2>REFERENCES</H2>
<P>Hofierka J., 1997. Direct solar radiation modelling within an
open GIS environment. Proceedings of JEC-GI'97 conference in Vienna,
@@ -155,6 +163,7 @@
Photovoltaic Assessments. <A HREF="http://www.blackwellpublishing.com/toc.asp?ref=1361-1682">Transactions
in GIS</A>, 8(2), 175-190</P>
+
<H2>AUTHORS</H2>
Thomas Huld, Joint Research Centre of
the European Commission, Ispra, Italy
@@ -174,4 +183,5 @@
<A HREF="mailto:Marcel.Suri at jrc.it">Marcel.Suri at jrc.it</A>
</ADDRESS>
-<P><I>Last changed: $Date: 2007/05/16 16:22:04 $</I>
+<P>
+<I>Last changed: $Date: 2007/05/16 16:22:04 $</I>
Modified: grass/branches/develbranch_6/raster/r.horizon/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.horizon/main.c 2009-08-11 11:50:41 UTC (rev 38684)
+++ grass/branches/develbranch_6/raster/r.horizon/main.c 2009-08-11 11:58:02 UTC (rev 38685)
@@ -49,7 +49,7 @@
#define SMALL 1.e-20
#define EPS 1.e-4
#define DIST "1.0"
-#define DEGREEINMETERS 111120.
+#define DEGREEINMETERS 111120. /* 1852m/nm * 60nm/degree = 111120 m/deg */
#define TANMINANGLE 0.008727 /* tan of minimum horizon angle (0.5 deg) */
#define AMAX1(arg1, arg2) ((arg1) >= (arg2) ? (arg1) : (arg2))
@@ -71,7 +71,7 @@
char *horizon = NULL;
char *mapset = NULL;
char *per;
-char shad_filename[256];
+char shad_filename[GNAME_MAX];
struct Cell_head cellhd;
struct Key_value *in_proj_info, *in_unit_info;
@@ -108,7 +108,7 @@
int ip, jp, ip100, jp100;
int n, m, m100, n100;
-int degreeOutput = 0;
+int degreeOutput = FALSE;
float **z, **z100, **horizon_raster;
double stepx, stepy, stepxhalf, stepyhalf, stepxy, xp, yp, op, dp, xg0, xx0,
yg0, yy0, deltx, delty;
@@ -139,9 +139,10 @@
mode = val;
}
-int ll_correction = 0;
+int ll_correction = FALSE;
double coslatsq;
+/* use G_distance() instead ??!?! */
double distance(double x1, double x2, double y1, double y2)
{
if (ll_correction) {
@@ -339,9 +340,8 @@
else {
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.");
-
+ G_fatal_error(
+ _("Can't read the coordinates from the \"coord\" option."));
}
/* Transform the coordinates to row/column */
@@ -358,15 +358,13 @@
if (isMode(WHOLE_RASTER)) {
if ((parm.direction->answer == NULL) && (parm.step->answer == NULL)) {
- G_fatal_error
- (_("You didn't specify a direction value or step size. Aborting."));
+ G_fatal_error(
+ _("You didn't specify a direction value or step size. Aborting."));
}
-
if (parm.horizon->answer == NULL) {
- G_fatal_error
- (_("You didn't specify a horizon raster name. Aborting."));
-
+ G_fatal_error(
+ _("You didn't specify a horizon raster name. Aborting."));
}
horizon = parm.horizon->answer;
if (parm.step->answer != NULL)
@@ -375,9 +373,8 @@
else {
if (parm.step->answer == NULL) {
- G_fatal_error
- (_("You didn't specify an angle step size. Aborting."));
-
+ G_fatal_error(
+ _("You didn't specify an angle step size. Aborting."));
}
sscanf(parm.step->answer, "%lf", &step);
@@ -389,43 +386,39 @@
}
if (parm.bufferzone->answer != NULL) {
- if (sscanf(parm.bufferzone->answer, "%lf", &bufferZone) != 1) {
+ if (sscanf(parm.bufferzone->answer, "%lf", &bufferZone) != 1)
G_fatal_error(_("Could not read bufferzone size. Aborting."));
- }
}
if (parm.e_buff->answer != NULL) {
- if (sscanf(parm.e_buff->answer, "%lf", &ebufferZone) != 1) {
- G_fatal_error(_("Could not read east bufferzone size. Aborting."));
- }
+ if (sscanf(parm.e_buff->answer, "%lf", &ebufferZone) != 1)
+ G_fatal_error(_("Could not read %s bufferzone size. Aborting."),
+ _("east"));
}
if (parm.w_buff->answer != NULL) {
- if (sscanf(parm.w_buff->answer, "%lf", &wbufferZone) != 1) {
- G_fatal_error(_("Could not read west bufferzone size. Aborting."));
- }
+ if (sscanf(parm.w_buff->answer, "%lf", &wbufferZone) != 1)
+ G_fatal_error(_("Could not read %s bufferzone size. Aborting."),
+ _("west"));
}
if (parm.s_buff->answer != NULL) {
- if (sscanf(parm.s_buff->answer, "%lf", &sbufferZone) != 1) {
- G_fatal_error
- (_("Could not read south bufferzone size. Aborting."));
- }
+ if (sscanf(parm.s_buff->answer, "%lf", &sbufferZone) != 1)
+ G_fatal_error(
+ _("Could not read %s bufferzone size. Aborting."),
+ _("south"));
}
-
-
if (parm.n_buff->answer != NULL) {
- if (sscanf(parm.n_buff->answer, "%lf", &nbufferZone) != 1) {
- G_fatal_error
- (_("Could not read north bufferzone size. Aborting."));
- }
+ if (sscanf(parm.n_buff->answer, "%lf", &nbufferZone) != 1)
+ G_fatal_error(
+ _("Could not read %s bufferzone size. Aborting."),
+ _("north"));
}
if (parm.maxdistance->answer != NULL) {
- if (sscanf(parm.maxdistance->answer, "%lf", &fixedMaxLength) != 1) {
+ if (sscanf(parm.maxdistance->answer, "%lf", &fixedMaxLength) != 1)
G_fatal_error(_("Could not read maximum distance. Aborting."));
- }
}
@@ -463,13 +456,17 @@
deltx = fabs(new_cellhd.east - new_cellhd.west);
delty = fabs(new_cellhd.north - new_cellhd.south);
- n /*n_cols */ = new_cellhd.cols;
- m /*n_rows */ = new_cellhd.rows;
- /*G_debug(3,"%lf %lf %lf %lf \n",ymax, ymin, xmin,xmax);
- */
+ n /* n_cols */ = new_cellhd.cols;
+ m /* n_rows */ = new_cellhd.rows;
+ /* G_debug(3,"%lf %lf %lf %lf \n",ymax, ymin, xmin,xmax); */
n100 = ceil(n / 100.);
m100 = ceil(m / 100.);
+/*
+ printf(". n=%f s=%f\n",new_cellhd.north, new_cellhd.south);
+ printf(". e=%f w=%f\n",new_cellhd.east, new_cellhd.west);
+ printf(". r=%d c=%d\n",new_cellhd.rows, new_cellhd.cols);
+*/
if (G_set_window(&new_cellhd) == -1)
exit(EXIT_FAILURE);
}
@@ -477,15 +474,16 @@
struct Key_Value *in_proj_info, *in_unit_info;
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!"));
+ G_fatal_error(
+ _("Can't get projection info of current location: "
+ "please set latitude via 'lat' or 'latin' option!"));
if ((in_unit_info = G_get_projunits()) == NULL)
G_fatal_error(_("Can't get projection units of current location"));
if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
- G_fatal_error
- (_("Can't get projection key values of current location"));
+ G_fatal_error(
+ _("Can't get projection key values of current location"));
G_free_key_value(in_proj_info);
G_free_key_value(in_unit_info);
@@ -498,8 +496,6 @@
G_fatal_error(_("Unable to set up lat/long projection parameters"));
-
-
/**********end of parser - ******************************/
@@ -515,7 +511,7 @@
exit(EXIT_FAILURE);
}
- /*sorry, I've moved OUTGR to calculate() - into the loop */
+ /* sorry, I've moved OUTGR() to calculate() - into the loop */
/* if(isMode(WHOLE_RASTER))
{
OUTGR(cellhd.rows,cellhd.cols);
@@ -529,6 +525,7 @@
/**********************end of main.c*****************/
+
int INPUT(void)
{
FCELL *cell1;
@@ -569,7 +566,7 @@
}
G_close_cell(fd1);
- /*create low resolution array 100 */
+ /* create low resolution array 100 */
for (i = 0; i < m100; i++) {
lmax = (i + 1) * 100;
if (lmax > m)
@@ -591,11 +588,10 @@
}
- /*find max Z */
+ /* find max Z */
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
zmax = amax1(zmax, z[i][j]);
-
}
}
@@ -619,10 +615,9 @@
cell1 = G_allocate_f_raster_buf();
fd1 = G_open_fp_cell_new(shad_filename);
if (fd1 < 0)
- G_fatal_error(_("Unable to create raster map %s"), shad_filename);
+ G_fatal_error(_("Unable to create raster map <%s>"), shad_filename);
}
-
if (numrows != G_window_rows())
G_fatal_error(_("OOPS: rows changed from %d to %d"), numrows,
G_window_rows());
@@ -643,13 +638,10 @@
}
G_put_f_raster_row(fd1, cell1);
}
-
} /* End loop over rows. */
G_close_cell(fd1);
-
-
return 1;
}
@@ -761,7 +753,6 @@
double calculate_shadow_onedirection(double shadow_angle)
{
-
shadow_angle = horizon_height();
return shadow_angle;
@@ -970,20 +961,20 @@
}
mindel = min(delx, dely);
- /*G_debug(3,"%d %d %d %lf %lf\n",ip, jp, mindel,xg0, yg0);*/
+ /* G_debug(3,"%d %d %d %lf %lf\n",ip, jp, mindel,xg0, yg0);*/
yy0 = yy0 + (mindel * stepsinangle);
xx0 = xx0 + (mindel * stepcosangle);
- /*G_debug(3," %lf %lf\n",xx0,yy0);*/
+ /* G_debug(3," %lf %lf\n",xx0,yy0);*/
return (3);
}
else {
- return (1); /*change of low res array - new cell is reaching limit for high resolution processing */
+ return (1); /* change of low res array - new cell is reaching limit for high resolution processing */
}
}
else {
- return (1); /*no change of low res array */
+ return (1); /* no change of low res array */
}
}
@@ -1003,9 +994,8 @@
if (succes != 1) {
break;
}
- /*
- curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS));
- */
+
+ /* curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS)); */
curvature_diff = 0.5 * length * length * invEarth;
z2 = z_orig + curvature_diff + length * tanh0;
@@ -1050,8 +1040,8 @@
double delt_dist;
char formatString[10];
+ char msg_buff[256];
-
int hor_row_start = buffer_s;
int hor_row_end = m - buffer_n;
@@ -1068,7 +1058,7 @@
yindex = (int)((ycoord - ymin) / stepy);
if ((G_projection() == PROJECTION_LL)) {
- ll_correction = 1;
+ ll_correction = TRUE;
}
@@ -1116,9 +1106,8 @@
horizon_raster[j][i] = 0.;
}
}
- /*
- definition of horizon angle in loop
- */
+
+ /* definition of horizon angle in loop */
if (step == 0.0) {
dfr_rad = 0;
arrayNumInt = 1;
@@ -1137,12 +1126,13 @@
if (step != 0.0)
sprintf(shad_filename, formatString, horizon, k);
+
angle = (single_direction * deg2rad) + (dfr_rad * k);
/*
com_par(angle);
*/
- G_message(_("Calculating map %01d of %01d (angle %lf, raster map <%s>)"), (k + 1), arrayNumInt,
- angle * rad2deg, shad_filename);
+ G_message(_("Calculating map %01d of %01d (angle %.2f, raster map <%s>)"),
+ (k + 1), arrayNumInt, angle * rad2deg, shad_filename);
for (j = hor_row_start; j < hor_row_end; j++) {
G_percent(j - hor_row_start, hor_numrows - 1, 2);
@@ -1168,11 +1158,10 @@
if ((G_projection() != PROJECTION_LL)) {
-
- if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0)
G_fatal_error("Error in pj_do_proj");
- }
}
+
latitude *= deg2rad;
longitude *= deg2rad;
@@ -1182,16 +1171,15 @@
twopi) ? inputAngle - twopi : inputAngle;
- delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
+ delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
latitude = (latitude + delt_lat) * rad2deg;
longitude = (longitude + delt_lon) * rad2deg;
if ((G_projection() != PROJECTION_LL)) {
- if (pj_do_proj(&longitude, &latitude, &oproj, &iproj) < 0) {
+ if (pj_do_proj(&longitude, &latitude, &oproj, &iproj) < 0)
G_fatal_error("Error in pj_do_proj");
- }
}
delt_east = longitude - xp;
@@ -1227,11 +1215,12 @@
maxlength =
(maxlength <
fixedMaxLength) ? maxlength : fixedMaxLength;
+
if (z_orig != UNDEFZ) {
- /*G_debug(3,"**************new line %d %d\n", i, j);
- */
+ /* G_debug(3,"**************new line %d %d\n", i, j); */
shadow_angle = horizon_height();
+
if (degreeOutput) {
shadow_angle *= rad2deg;
}
@@ -1246,7 +1235,7 @@
horizon_raster[j - buffer_s][i - buffer_w] =
shadow_angle;
- } /* undefs */
+ } /* undefs */
}
}
@@ -1258,19 +1247,38 @@
horizon_raster[j][i] = 0.;
}
- /*return back the buffered region */
+ /* return back the buffered region */
if (bufferZone > 0.) {
if (G_set_window(&new_cellhd) == -1)
exit(0);
}
+ /* write metadata */
G_short_history(shad_filename, "raster", &history);
+
+ sprintf(msg_buff,
+ "Angular height of terrain horizon, map %01d of %01d",
+ (k + 1), arrayNumInt);
+ G_put_cell_title(shad_filename, msg_buff);
+
+ if (degreeOutput)
+ G_write_raster_units(shad_filename, "degrees");
+ else
+ G_write_raster_units(shad_filename, "radians");
+
G_command_history(&history);
- G_write_history(shad_filename, &history);
+ /* insert a blank line */
+ history.edhist[history.edlinecnt][0] = '\0';
+ history.edlinecnt++;
+ sprintf(msg_buff,
+ "Horizon view from azimuth angle %.2f degrees CCW from East",
+ angle * rad2deg);
+ strcpy(history.edhist[history.edlinecnt], msg_buff);
+ history.edlinecnt++;
+
+ G_write_history(shad_filename, &history);
}
-
}
-
}
More information about the grass-commit
mailing list