[GRASS-SVN] r38686 - grass/trunk/raster/r.horizon

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Aug 11 08:06:15 EDT 2009


Author: hamish
Date: 2009-08-11 08:06:14 -0400 (Tue, 11 Aug 2009)
New Revision: 38686

Modified:
   grass/trunk/raster/r.horizon/main.c
   grass/trunk/raster/r.horizon/r.horizon.html
Log:
write out metadata; whitespace/msg cleanup (merge from devbr6)

Modified: grass/trunk/raster/r.horizon/main.c
===================================================================
--- grass/trunk/raster/r.horizon/main.c	2009-08-11 11:58:02 UTC (rev 38685)
+++ grass/trunk/raster/r.horizon/main.c	2009-08-11 12:06:14 UTC (rev 38686)
@@ -50,7 +50,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))
@@ -72,7 +72,7 @@
 const char *horizon = NULL;
 const char *mapset = NULL;
 const char *per;
-char shad_filename[256];
+char shad_filename[GNAME_MAX];
 
 struct Cell_head cellhd;
 struct Key_value *in_proj_info, *in_unit_info;
@@ -109,7 +109,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;
@@ -140,9 +140,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) {
@@ -340,9 +341,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 */
@@ -359,15 +359,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)
@@ -376,9 +374,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);
 
@@ -390,43 +387,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."));
-	}
     }
 
 
@@ -464,10 +457,9 @@
 	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.);
 
@@ -478,15 +470,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);
@@ -499,8 +492,6 @@
 	G_fatal_error(_("Unable to set up lat/long projection parameters"));
 
 
-
-
 /**********end of parser - ******************************/
 
 
@@ -516,7 +507,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);
@@ -530,6 +521,7 @@
 
 /**********************end of main.c*****************/
 
+
 int INPUT(void)
 {
     FCELL *cell1;
@@ -570,7 +562,7 @@
     }
     Rast_close(fd1);
 
-    /*create low resolution array 100 */
+    /* create low resolution array 100 */
     for (i = 0; i < m100; i++) {
 	lmax = (i + 1) * 100;
 	if (lmax > m)
@@ -592,11 +584,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]);
-
 	}
     }
 
@@ -620,10 +611,9 @@
 	cell1 = Rast_allocate_f_buf();
 	fd1 = Rast_open_fp_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());
@@ -644,13 +634,10 @@
 	    }
 	    Rast_put_f_row(fd1, cell1);
 	}
-
     }				/* End loop over rows. */
 
     Rast_close(fd1);
 
-
-
     return 1;
 }
 
@@ -762,7 +749,6 @@
 
 double calculate_shadow_onedirection(double shadow_angle)
 {
-
     shadow_angle = horizon_height();
 
     return shadow_angle;
@@ -971,20 +957,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 */
     }
 }
 
@@ -1004,9 +990,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;
@@ -1051,8 +1036,8 @@
     double delt_dist;
 
     char formatString[10];
+    char msg_buff[256];
 
-
     int hor_row_start = buffer_s;
     int hor_row_end = m - buffer_n;
 
@@ -1069,7 +1054,7 @@
     yindex = (int)((ycoord - ymin) / stepy);
 
     if ((G_projection() == PROJECTION_LL)) {
-	ll_correction = 1;
+	ll_correction = TRUE;
     }
 
 
@@ -1117,9 +1102,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;
@@ -1138,12 +1122,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);
@@ -1169,11 +1154,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;
 
@@ -1183,16 +1167,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;
@@ -1228,11 +1211,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;
 			}
@@ -1247,7 +1231,7 @@
 			horizon_raster[j - buffer_s][i - buffer_w] =
 			    shadow_angle;
 
-		    }		/* undefs */
+		    }	/* undefs */
 		}
 	    }
 
@@ -1259,19 +1243,38 @@
 		    horizon_raster[j][i] = 0.;
 	    }
 
-	    /*return back the buffered region */
+	    /* return back the buffered region */
 	    if (bufferZone > 0.) {
 		if (Rast_set_window(&new_cellhd) == -1)
 		    exit(0);
 	    }
 
+	    /* write metadata */
 	    Rast_short_history(shad_filename, "raster", &history);
+
+	    sprintf(msg_buff,
+		    "Angular height of terrain horizon, map %01d of %01d",
+		    (k + 1), arrayNumInt);
+	    Rast_put_cell_title(shad_filename, msg_buff);
+
+	    if (degreeOutput)
+		Rast_write_units(shad_filename, "degrees");
+	    else
+		Rast_write_units(shad_filename, "radians");
+
 	    Rast_command_history(&history);
-	    Rast_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++;
+
+	    Rast_write_history(shad_filename, &history);
 	}
-
     }
-
 }

Modified: grass/trunk/raster/r.horizon/r.horizon.html
===================================================================
--- grass/trunk/raster/r.horizon/r.horizon.html	2009-08-11 11:58:02 UTC (rev 38685)
+++ grass/trunk/raster/r.horizon/r.horizon.html	2009-08-11 12:06:14 UTC (rev 38686)
@@ -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)</dd>
 </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 &gt; horizon.out
+r.horizon elevin=DEM horizonstep=30 direction=15 bufferzone=200 \
+    coord=47.302,7.365 dist=0.7 &gt; 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>
 <P>
 Thomas Huld, Joint Research Centre of
@@ -178,4 +187,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> 



More information about the grass-commit mailing list