[GRASS-SVN] r51947 - in grass/trunk/raster: . r.shaded.relief2
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Jun 3 04:29:43 PDT 2012
Author: mmetz
Date: 2012-06-03 04:29:43 -0700 (Sun, 03 Jun 2012)
New Revision: 51947
Added:
grass/trunk/raster/r.shaded.relief2/
grass/trunk/raster/r.shaded.relief2/main.c
grass/trunk/raster/r.shaded.relief2/r.shaded.relief2.html
Removed:
grass/trunk/raster/r.shaded.relief2/r.shaded.relief.html
grass/trunk/raster/r.shaded.relief2/r.shaded.relief.py
Modified:
grass/trunk/raster/r.shaded.relief2/Makefile
Log:
r.shaded.relief2: enhanced fast C version of r.shaded.relief
Modified: grass/trunk/raster/r.shaded.relief2/Makefile
===================================================================
--- grass/trunk/scripts/r.shaded.relief/Makefile 2012-06-02 17:12:24 UTC (rev 51942)
+++ grass/trunk/raster/r.shaded.relief2/Makefile 2012-06-03 11:29:43 UTC (rev 51947)
@@ -1,7 +1,10 @@
MODULE_TOPDIR = ../..
-PGM=r.shaded.relief
+PGM = r.shaded.relief2
-include $(MODULE_TOPDIR)/include/Make/Script.make
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
-default: script
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass/trunk/raster/r.shaded.relief2/main.c
===================================================================
--- grass/trunk/raster/r.shaded.relief2/main.c (rev 0)
+++ grass/trunk/raster/r.shaded.relief2/main.c 2012-06-03 11:29:43 UTC (rev 51947)
@@ -0,0 +1,463 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.shaded.relief
+ * AUTHOR(S): CERL
+ * parameters standardized: Markus Neteler, 2008
+ * updates: Michael Barton, 2004
+ * updates: Gordon Keith, 2003
+ * updates: Andreas Lange, 2001
+ * updates: David Finlayson, 2001
+ * updates: Markus Neteler, 2001, 1999
+ * Converted to Python by Glynn Clements
+ * Converted to C by Markus Metz
+ * PURPOSE: Creates shaded relief map from raster elevation map (DEM)
+ * COPYRIGHT: (C) 1999 - 2008, 2010, 2012 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+/*
+ * July 2007 - allow input from other mapsets (Brad Douglas)
+ *
+ * May 2005 - fixed wrong units parameter (Markus Neteler)
+ *
+ * September 2004 - Added z exaggeration control (Michael Barton)
+ * April 2004 - updated for GRASS 5.7 by Michael Barton
+ *
+ * 9/2004 Adds scale factor input (as per documentation); units set scale only if specified for lat/long regions
+ * Also, adds option of controlling z-exaggeration.
+ *
+ * 6/2003 fixes for Lat/Long Gordon Keith <gordon.keith at csiro.au>
+ * If n is a number then the ewres and nsres are mulitplied by that scale
+ * to calculate the shading.
+ * If n is the letter M (either case) the number of metres is degree of
+ * latitude is used as the scale.
+ * If n is the letter f then the number of feet in a degree is used.
+ * It scales latitude and longitude equally, so it's only approximately
+ * right, but for shading its close enough. It makes the difference
+ * between an unusable and usable shade.
+ *
+ * 10/2001 fix for testing for dashes in raster file name
+ * by Andreas Lange <andreas.lange at rhein-main.de>
+ * 10/2001 added parser support - Markus Neteler
+ * 9/2001 fix to keep NULLs as is (was value 22 before) - Markus Neteler
+ * 1/2001 fix for NULL by David Finlayson <david_finlayson at yahoo.com>
+ * 11/99 updated $ewres to ewres() and $nsres to nsres()
+ * updated number to FP in r.mapcalc statement Markus Neteler
+ */
+
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+int main(int argc, char *argv[])
+{
+ int in_fd;
+ int out_fd;
+ DCELL *elev_cell[3], *temp;
+ DCELL *out_rast, *out_ptr = NULL;
+ DCELL *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9;
+ int Wrap; /* global wraparound */
+ struct Cell_head window;
+ struct History hist;
+ struct Colors colors;
+
+ const char *elev_name;
+ const char *sr_name;
+ int out_type = CELL_TYPE;
+ size_t out_size;
+ const char *units;
+ char buf[GNAME_MAX];
+ int nrows, row;
+ int ncols, col;
+
+ double zmult, scale, altitude, azimuth;
+ double north, east, south, west, ns_med;
+
+ double degrees_to_radians, radians_to_degrees;
+ double H, V;
+ double dx; /* partial derivative in ew direction */
+ double dy; /* partial derivative in ns direction */
+ double key;
+ double slp_in_rad, aspect, cang;
+
+ struct FPRange range;
+ DCELL min, max;
+
+ struct GModule *module;
+ struct
+ {
+ struct Option *elevation, *relief, *altitude, *azimuth, *zmult,
+ *scale, *units;
+ } parm;
+ char *desc;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("terrain"));
+ module->label = _("Creates shaded relief map from an elevation map (DEM).");
+
+ parm.elevation = G_define_standard_option(G_OPT_R_INPUT);
+
+ parm.relief = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.relief->required = NO;
+ parm.relief->label = _("Name for output shaded relief map");
+ parm.relief->description = _("Default: <input_map>.shade");
+
+ parm.altitude = G_define_option();
+ parm.altitude->key = "altitude";
+ parm.altitude->type = TYPE_DOUBLE;
+ parm.altitude->required = NO;
+ parm.altitude->answer = "30";
+ parm.altitude->options = "0-90";
+ parm.altitude->description = _("Altitude of the sun in degrees above the horizon");
+
+ parm.azimuth = G_define_option();
+ parm.azimuth->key = "azimuth";
+ parm.azimuth->type = TYPE_DOUBLE;
+ parm.azimuth->required = NO;
+ parm.azimuth->answer = "270";
+ parm.azimuth->options = "0-360";
+ parm.azimuth->description =
+ _("Azimuth of the sun in degrees to the east of north");
+
+ parm.zmult = G_define_option();
+ parm.zmult->key = "zmult";
+ parm.zmult->type = TYPE_DOUBLE;
+ parm.zmult->required = NO;
+ parm.zmult->answer = "1";
+ parm.zmult->description = _("Factor for exaggerating relief");
+
+ parm.scale = G_define_option();
+ parm.scale->key = "scale";
+ parm.scale->type = TYPE_DOUBLE;
+ parm.scale->required = NO;
+ parm.scale->answer = "1";
+ parm.scale->description =
+ _("Scale factor for converting meters to elevation units");
+
+ parm.units = G_define_option();
+ parm.units->key = "units";
+ parm.units->type = TYPE_STRING;
+ parm.units->required = NO;
+ parm.units->options = "intl,survey";
+ parm.units->description = _("Elevation units (overrides scale factor)");
+ desc = NULL;
+ G_asprintf(&desc,
+ "intl;%s;"
+ "survey;%s",
+ _("international feet"),
+ _("survey feet"));
+ parm.units->descriptions = desc;
+
+
+ degrees_to_radians = M_PI / 180.0;
+ radians_to_degrees = 180. / M_PI;
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ elev_name = parm.elevation->answer;
+ if (parm.relief->answer) {
+ sr_name = parm.relief->answer;
+ }
+ else {
+ char xname[GNAME_MAX], xmapset[GNAME_MAX];
+
+ if (G_name_is_fully_qualified(elev_name, xname, xmapset))
+ sprintf(buf, "%s.shade", xname);
+ else
+ sprintf(buf, "%s.shade", elev_name);
+ sr_name = G_store(buf);
+ }
+
+ G_check_input_output_name(elev_name, sr_name, G_FATAL_EXIT);
+
+ if (sscanf(parm.altitude->answer, "%lf", &altitude) != 1 || altitude < 0.0) {
+ G_fatal_error(_("%s=%s - must be a non-negative number"),
+ parm.altitude->key, parm.altitude->answer);
+ }
+ altitude *= degrees_to_radians;
+
+ if (sscanf(parm.azimuth->answer, "%lf", &azimuth) != 1 || azimuth < 0.0) {
+ G_fatal_error(_("%s=%s - must be a non-negative number"),
+ parm.azimuth->key, parm.azimuth->answer);
+ }
+ /* correct azimuth to East (GRASS convention):
+ * this seems to be backwards, but in fact it works so leave it. */
+ azimuth = (azimuth - 90.) * degrees_to_radians;
+
+ if (sscanf(parm.zmult->answer, "%lf", &zmult) != 1 || zmult == 0.0) {
+ G_fatal_error(_("%s=%s - must not be zero"),
+ parm.zmult->key, parm.zmult->answer);
+ }
+
+ if (sscanf(parm.scale->answer, "%lf", &scale) != 1 || scale <= 0.0) {
+ G_fatal_error(_("%s=%s - must be a positive number"),
+ parm.scale->key,
+ parm.scale->answer);
+ }
+
+ G_get_set_window(&window);
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* horizontal distances are calculated in meters by G_distance() */
+ if (parm.units->answer) {
+ units = parm.units->answer;
+ if (strcmp(units, "intl") == 0) {
+ /* 1 international foot = 0.3048 meters */
+ scale = 1. / 0.3048;
+
+ }
+ else if (strcmp(units, "survey") == 0) {
+ /* 1 survey foot = 1200 / 3937 meters */
+ scale = 3937. / 1200.;
+ }
+ }
+
+ Wrap = 0;
+ if (G_projection() == PROJECTION_LL) {
+ if ((window.west == (window.east - 360.))
+ || (window.east == (window.west - 360.))) {
+ Wrap = 1;
+ ncols += 2;
+ }
+ }
+
+ /* H = window.ew_res * 4 * 2/ zmult; *//* horizontal (east-west) run
+ times 4 for weighted difference */
+ /* V = window.ns_res * 4 * 2/ zmult; *//* vertical (north-south) run
+ times 4 for weighted difference */
+
+
+ G_begin_distance_calculations();
+ north = Rast_row_to_northing(0.5, &window);
+ ns_med = Rast_row_to_northing(1.5, &window);
+ south = Rast_row_to_northing(2.5, &window);
+ east = Rast_col_to_easting(2.5, &window);
+ west = Rast_col_to_easting(0.5, &window);
+ V = G_distance(east, north, east, south) * 4 * scale / zmult;
+ H = G_distance(east, ns_med, west, ns_med) * 4 * scale / zmult;
+ /* ____________________________
+ |c1 |c2 |c3 |
+ | | | |
+ | | north | |
+ | | | |
+ |________|________|________|
+ |c4 |c5 |c6 |
+ | | | |
+ | east | ns_med | west |
+ | | | |
+ |________|________|________|
+ |c7 |c8 |c9 |
+ | | | |
+ | | south | |
+ | | | |
+ |________|________|________|
+ */
+
+ /* open the elevation file for reading */
+ in_fd = Rast_open_old(elev_name, "");
+ elev_cell[0] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
+ Rast_set_d_null_value(elev_cell[0], ncols);
+ elev_cell[1] = (DCELL *) G_calloc(ncols, sizeof(DCELL));
+ Rast_set_d_null_value(elev_cell[1], ncols);
+ elev_cell[2] = (DCELL *) G_calloc(ncols, sizeof(DCELL));
+ Rast_set_d_null_value(elev_cell[2], ncols);
+
+ out_fd = Rast_open_new(sr_name, out_type);
+ out_rast = Rast_allocate_buf(out_type);
+ Rast_set_null_value(out_rast, Rast_window_cols(), out_type);
+ Rast_put_row(out_fd, out_rast, out_type);
+ out_size = Rast_cell_size(out_type);
+
+ if (Wrap) {
+ Rast_get_d_row_nomask(in_fd, elev_cell[1] + 1, 0);
+ elev_cell[1][0] = elev_cell[1][Rast_window_cols() - 1];
+ elev_cell[1][Rast_window_cols() + 1] = elev_cell[1][2];
+ }
+ else
+ Rast_get_d_row_nomask(in_fd, elev_cell[1], 0);
+
+ if (Wrap) {
+ Rast_get_d_row_nomask(in_fd, elev_cell[2] + 1, 1);
+ elev_cell[2][0] = elev_cell[2][Rast_window_cols() - 1];
+ elev_cell[2][Rast_window_cols() + 1] = elev_cell[2][2];
+ }
+ else
+ Rast_get_d_row_nomask(in_fd, elev_cell[2], 1);
+
+ G_verbose_message(_("Percent complete..."));
+
+ for (row = 2; row < nrows; row++) {
+ /* if projection is Lat/Lon, recalculate V and H */
+ if (G_projection() == PROJECTION_LL) {
+ north = Rast_row_to_northing((row - 2 + 0.5), &window);
+ ns_med = Rast_row_to_northing((row - 1 + 0.5), &window);
+ south = Rast_row_to_northing((row + 0.5), &window);
+ east = Rast_col_to_easting(2.5, &window);
+ west = Rast_col_to_easting(0.5, &window);
+ V = G_distance(east, north, east, south) * 4 * scale / zmult;
+ H = G_distance(east, ns_med, west, ns_med) * 4 * scale / zmult;
+ /* ____________________________
+ |c1 |c2 |c3 |
+ | | | |
+ | | north | |
+ | | | |
+ |________|________|________|
+ |c4 |c5 |c6 |
+ | | | |
+ | east | ns_med | west |
+ | | | |
+ |________|________|________|
+ |c7 |c8 |c9 |
+ | | | |
+ | | south | |
+ | | | |
+ |________|________|________|
+ */
+ }
+
+ G_percent(row, nrows, 2);
+ temp = elev_cell[0];
+ elev_cell[0] = elev_cell[1];
+ elev_cell[1] = elev_cell[2];
+ elev_cell[2] = temp;
+
+ if (Wrap) {
+ Rast_get_d_row_nomask(in_fd, elev_cell[2] + 1, row);
+ elev_cell[2][0] = elev_cell[2][Rast_window_cols() - 1];
+ elev_cell[2][Rast_window_cols() + 1] = elev_cell[2][2];
+ }
+ else
+ Rast_get_d_row_nomask(in_fd, elev_cell[2], row);
+
+ c1 = elev_cell[0];
+ c2 = c1 + 1;
+ c3 = c1 + 2;
+ c4 = elev_cell[1];
+ c5 = c4 + 1;
+ c6 = c4 + 2;
+ c7 = elev_cell[2];
+ c8 = c7 + 1;
+ c9 = c7 + 2;
+
+ if (Wrap)
+ out_ptr = out_rast;
+ else
+ out_ptr = G_incr_void_ptr(out_rast, out_size);
+
+ /*skip first cell of the row */
+
+ for (col = ncols - 2; col-- > 0;
+ c1++, c2++, c3++, c4++, c5++, c6++, c7++, c8++, c9++) {
+ /* DEBUG:
+ fprintf(stdout, "\n%.0f %.0f %.0f\n%.0f %.0f %.0f\n%.0f %.0f %.0f\n",
+ *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9);
+ */
+
+ if (Rast_is_d_null_value(c1) || Rast_is_d_null_value(c2) ||
+ Rast_is_d_null_value(c3) || Rast_is_d_null_value(c4) ||
+ Rast_is_d_null_value(c5) || Rast_is_d_null_value(c6) ||
+ Rast_is_d_null_value(c7) || Rast_is_d_null_value(c8) ||
+ Rast_is_d_null_value(c9)) {
+
+ Rast_set_null_value(out_ptr, 1, out_type);
+ out_ptr = G_incr_void_ptr(out_ptr, out_size);
+
+ continue;
+ } /* no data */
+
+ /* shaded relief */
+ /* slope */
+ dx = (*c1 + 2 * *c4 + *c7 - *c3 - 2 * *c6 - *c9) / H;
+ dy = (*c1 + 2 * *c2 + *c3 - *c7 - 2 * *c8 - *c9) / V;
+
+ key = dx * dx + dy * dy;
+
+ slp_in_rad = M_PI / 2. - atan(sqrt(key));
+
+ /* aspect */
+ aspect = atan2(dy, dx);
+
+ if (aspect != aspect)
+ aspect = degrees_to_radians;
+ if (dx != 0 || dy != 0) {
+ if (aspect == 0)
+ aspect = 2 * M_PI;
+ }
+
+#if 0
+ /* the original script was rounding aspect. Why? */
+ aspect *= radians_to_degrees;
+ if (aspect < 0)
+ aspect = (int)(aspect - 0.5);
+ else
+ aspect = (int)(aspect + 0.5);
+ aspect *= degrees_to_radians;
+#endif
+
+ /* shaded relief */
+ cang = sin(altitude) * sin(slp_in_rad) +
+ cos(altitude) * cos(slp_in_rad) * cos(azimuth - aspect);
+
+ Rast_set_d_value(out_ptr, (DCELL) 255 * cang, out_type);
+
+ out_ptr = G_incr_void_ptr(out_ptr, out_size);
+
+ } /* column for loop */
+
+ Rast_put_row(out_fd, out_rast, out_type);
+
+ } /* row loop */
+
+ G_percent(row, nrows, 2);
+
+ Rast_close(in_fd);
+
+ Rast_set_null_value(out_rast, Rast_window_cols(), out_type);
+ Rast_put_row(out_fd, out_rast, out_type);
+ Rast_close(out_fd);
+
+ G_debug(1, "Creating support files...");
+
+ /* write colors for shaded relief */
+ Rast_init_colors(&colors);
+ Rast_read_fp_range(sr_name, G_mapset(), &range);
+ Rast_get_fp_range_min_max(&range, &min, &max);
+ min -= 0.01;
+ max += 0.01;
+ Rast_make_grey_scale_fp_colors(&colors, min, max);
+ Rast_write_colors(sr_name, G_mapset(), &colors);
+
+ sprintf(buf, "Shaded relief of \"%s\"", elev_name);
+ Rast_put_cell_title(sr_name, buf);
+
+ /* writing history file */
+ Rast_short_history(sr_name, "raster", &hist);
+ Rast_append_format_history(&hist, "r.shaded.relief settings:");
+ Rast_append_format_history(&hist,
+ "altitude=%f azimuth=%f zmult=%f scale=%f",
+ altitude * radians_to_degrees,
+ azimuth * radians_to_degrees,
+ zmult, scale);
+ Rast_format_history(&hist, HIST_DATSRC_1, "raster elevation file %s", elev_name);
+ Rast_command_history(&hist);
+ Rast_write_history(sr_name, &hist);
+
+ G_message(_("Shaded relief raster map <%s> complete"), sr_name);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/raster/r.shaded.relief2/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Deleted: grass/trunk/raster/r.shaded.relief2/r.shaded.relief.html
===================================================================
--- grass/trunk/scripts/r.shaded.relief/r.shaded.relief.html 2012-06-02 17:12:24 UTC (rev 51942)
+++ grass/trunk/raster/r.shaded.relief2/r.shaded.relief.html 2012-06-03 11:29:43 UTC (rev 51947)
@@ -1,95 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>r.shaded.relief</em> creates a raster shaded relief map based on current
-resolution settings and on sun altitude, azimuth, and z-exaggeration values
-entered by the user. If no output shademap name is given, the new shaded
-relief map is named <em><input_map>.shade</em>.
-The map is assigned a grey-scale color table.
-
-<p>The module then prompts the user to enter values for:
-<ol>
-<li>The <b>altitude</b> of the sun in degrees above the horizon
-(a value between 0 and 90 degrees), and
-
-<li>The <b>azimuth</b> of the sun in degrees to the east of north
-(a value between 0 and 360 degrees).
-
-<li>The name of a raster map layer to provide elevation values for the
-shaded relief map. Typically, this would be a map layer of elevation;
-however, any raster map layer can be named.
-
-<li>The scaling parameter, which compensates for a different horizontal
-<b>scale</b> than vertical scale. If 'scale' is a number then the ewres
-and nsres are multiplied by that scale to calculate the shading. (Default=1.0
-for equivalent horizontal and vertical scales.)
-
-<li>For the special case when a latitude-longitude projection is used with an
-elevation map measured in meters (e.g., SRTM, ETOPO2 etc.) or feet, the
-<b>units</b> can be set to automatically set the horizontal scale to the the number of
-meters (scale=111120) or feet (scale=370400) in a degree of latitude. The script scales
-latitude and longitude equally, so it's only approximately right, but for shading
-it's close enough. It makes the difference between a usable and unusable shade.
-The <b>units</b> parameter overrides the <b>scale</b> parameter.
-
-<li>The <b>zmult</b> exaggeration factor that changes the apparent relief
-for the shaded relief map. This can be any positive (or negative) floating
-point value. (Default=1.0)
-</ol>
-
-Specifically, <em>r.shaded.relief</em> executes a <em>r.mapcalc</em>
-statement. Refer to the manual entry for <em>r.mapcalc</em> for an explanation
-of the filtering syntax shown in the above expression.
-See, for example, the section on "The Neighborhood Modifier".
-
-<p><em>r.shaded.relief</em> then runs <em>r.colors</em>
-to assign a grey-scale color table to the new shaded relief map.
-
-<h2>NOTES</h2>
-
-To visually improve the result of shade maps from low resolution elevation
-models, use <em>r.resamp.interp</em> with bilinear or bicubic method to
-resample the DEM at higher resolution. <em>r.shaded.relief</em> is then
-run on the resampled DEM.
-
-<h2>EXAMPLES</h2>
-
-In this example, the aspect map in the North Carolina sample
-dataset location is used to hillshade the elevation map:
-
-<div class="code"><pre>
-g.region rast=elevation -p
-r.shaded.relief input=elevation output=elevation.shaded
-</pre></div>
-
-<p>In Latitude-Longitude locations (or other non-metric locations), the
-<em>scale</em> factor has to be used:
-
-<div class="code"><pre>
-# Latitude-Longitude example
-r.shaded.relief input=srtm output=srtm.shaded scale=111120
-</pre></div>
-
-<h2>SEE ALSO</h2>
-
-"<i>r.mapcalc: An Algebra for GIS and Image Processing</i>", by Michael Shapiro
-and Jim Westervelt, U.S. Army Construction Engineering Research Laboratory
-(March/1991) (available from the GRASS web site).
-
-<p><!-- RGB version not ported to GRASS 6:
- <em><a href="shade.clr.sh.html">shade.clr.sh</a></em><br>
- -->
-<em>
-<a href="d.his.html">d.his</a>,
-<a href="g.region.html">g.region</a>,
-<a href="r.blend.html">r.blend</a>,
-<a href="r.colors.html">r.colors</a>,
-<a href="r.mapcalc.html">r.mapcalc</a>,
-<a href="r.resamp.interp.html">r.resamp.interp</a>
-</em>
-
-<h2>AUTHOR</h2>
-
-Jim Westervelt, U.S. Army Construction Engineering
-Research Laboratory
-
-<p><i>Last changed: $Date$</i>
Deleted: grass/trunk/raster/r.shaded.relief2/r.shaded.relief.py
===================================================================
--- grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py 2012-06-02 17:12:24 UTC (rev 51942)
+++ grass/trunk/raster/r.shaded.relief2/r.shaded.relief.py 2012-06-03 11:29:43 UTC (rev 51947)
@@ -1,179 +0,0 @@
-#!/usr/bin/env python
-
-############################################################################
-#
-# MODULE: r.shaded.relief
-# AUTHOR(S): CERL
-# parameters standardized: Markus Neteler, 2008
-# updates: Michael Barton, 2004
-# updates: Gordon Keith, 2003
-# updates: Andreas Lange, 2001
-# updates: David Finlayson, 2001
-# updates: Markus Neteler, 2001, 1999
-# Converted to Python by Glynn Clements
-# PURPOSE: Creates shaded relief map from raster elevation map (DEM)
-# COPYRIGHT: (C) 1999 - 2008, 2010 by the GRASS Development Team
-#
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-#
-# July 2007 - allow input from other mapsets (Brad Douglas)
-#
-# May 2005 - fixed wrong units parameter (Markus Neteler)
-#
-# September 2004 - Added z exaggeration control (Michael Barton)
-# April 2004 - updated for GRASS 5.7 by Michael Barton
-#
-# 9/2004 Adds scale factor input (as per documentation); units set scale only if specified for lat/long regions
-# Also, adds option of controlling z-exaggeration.
-#
-# 6/2003 fixes for Lat/Long Gordon Keith <gordon.keith at csiro.au>
-# If n is a number then the ewres and nsres are mulitplied by that scale
-# to calculate the shading.
-# If n is the letter M (either case) the number of metres is degree of
-# latitude is used as the scale.
-# If n is the letter f then the number of feet in a degree is used.
-# It scales latitude and longitude equally, so it's only approximately
-# right, but for shading its close enough. It makes the difference
-# between an unusable and usable shade.
-#
-# 10/2001 fix for testing for dashes in raster file name
-# by Andreas Lange <andreas.lange at rhein-main.de>
-# 10/2001 added parser support - Markus Neteler
-# 9/2001 fix to keep NULLs as is (was value 22 before) - Markus Neteler
-# 1/2001 fix for NULL by David Finlayson <david_finlayson at yahoo.com>
-# 11/99 updated $ewres to ewres() and $nsres to nsres()
-# updated number to FP in r.mapcalc statement Markus Neteler
-
-#%module
-#% description: Creates shaded relief map from an elevation map (DEM).
-#% keywords: raster
-#% keywords: elevation
-#%end
-#%option G_OPT_R_INPUT
-#% description: Name of input elevation raster map
-#%end
-#%option G_OPT_R_OUTPUT
-#%end
-#%option
-#% key: altitude
-#% type: double
-#% description: Altitude of the sun in degrees above the horizon
-#% required: no
-#% options : 0-90
-#% answer: 30
-#%end
-#%option
-#% key: azimuth
-#% type: double
-#% description: Azimuth of the sun in degrees to the east of north
-#% required: no
-#% options : 0-360
-#% answer: 270
-#%end
-#%option
-#% key: zmult
-#% type: double
-#% description: Factor for exaggerating relief
-#% required: no
-#% answer: 1
-#%end
-#%option
-#% key: scale
-#% type: double
-#% description: Scale factor for converting horizontal units to elevation units
-#% required: no
-#% answer: 1
-#% guisection: Scaling
-#%end
-#%option
-#% key: units
-#% type: string
-#% description: Set scaling factor (applies to lat./long. locations only, none: scale=1)
-#% required: no
-#% options: none,meters,feet
-#% answer: none
-#% guisection: Scaling
-#%end
-
-import sys
-import os
-import string
-import grass.script as grass
-
-def main():
- input = options['input']
- output = options['output']
- altitude = options['altitude']
- azimuth = options['azimuth']
- zmult = options['zmult']
- scale = float(options['scale'])
- units = options['units']
- verbose_level = os.getenv('GRASS_VERBOSE')
- if verbose_level is None:
- verbose_level = 2
-
- if not grass.find_file(input)['file']:
- grass.fatal(_("Raster map <%s> not found") % input)
-
- if input == output:
- grass.fatal(_("Input elevation map and output relief map must have different names"))
-
- # LatLong locations only:
- if units == 'meters':
- # scale=111120
- scale *= 1852 * 60
-
- # LatLong locations only:
- if units == 'feet':
- # scale=364567.2
- scale *= 6076.12 * 60
-
- #correct azimuth to East (GRASS convention):
- # this seems to be backwards, but in fact it works so leave it.
- az = float(azimuth) - 90
-
- t = string.Template(
- r'''eval( \
- x=($zmult*$input[-1,-1] + 2*$zmult*$input[0,-1] + $zmult*$input[1,-1] - \
- $zmult*$input[-1,1] - 2*$zmult*$input[0,1] - $zmult*$input[1,1])/(8.*ewres()*$scale), \
- y=($zmult*$input[-1,-1] + 2*$zmult*$input[-1,0] + $zmult*$input[-1,1] - \
- $zmult*$input[1,-1] - 2*$zmult*$input[1,0] - $zmult*$input[1,1])/(8.*nsres()*$scale), \
- slope=90.-atan(sqrt(x*x + y*y)), \
- a=round(atan(x,y)), \
- a=if(isnull(a),1,a), \
- aspect=if(x!=0||y!=0,if(a,a,360.)), \
- cang = sin($altitude)*sin(slope) + cos($altitude)*cos(slope) * cos($az-aspect) \
- )
- $output = if(isnull(cang), null(), 100.*cang)''')
- expr = t.substitute(altitude = altitude, az = az, input = input, output = output, scale = scale, zmult = zmult)
- p = grass.feed_command('r.mapcalc')
- p.stdin.write(expr)
- p.stdin.close()
-
- if p.wait() != 0:
- grass.fatal(_("In calculation, script aborted."))
-
- if verbose_level < 3:
- grass.run_command('r.colors', map = output, color = 'grey', quiet = True)
- else:
- grass.run_command('r.colors', map = output, color = 'grey')
-
- # record metadata
- grass.run_command('r.support', map = output, title = 'Shaded relief of "%s"' % input, history = '')
- grass.run_command('r.support', map = output, history = "r.shaded.relief settings:")
- t = string.Template("altitude=$altitude azimuth=$azimuth zmult=$zmult scale=$scale")
- parms = dict(altitude = altitude, azimuth = azimuth, zmult = zmult, scale = options['scale'])
- grass.run_command('r.support', map = output, history = t.substitute(parms))
- if units:
- grass.run_command('r.support', map = output, history = " units=%s (adjusted scale=%s)" % (units, scale))
-
- # write cmd history:
- grass.raster_history(output)
-
-if __name__ == "__main__":
- options, flags = grass.parser()
- main()
Copied: grass/trunk/raster/r.shaded.relief2/r.shaded.relief2.html (from rev 51942, grass/trunk/scripts/r.shaded.relief/r.shaded.relief.html)
===================================================================
--- grass/trunk/raster/r.shaded.relief2/r.shaded.relief2.html (rev 0)
+++ grass/trunk/raster/r.shaded.relief2/r.shaded.relief2.html 2012-06-03 11:29:43 UTC (rev 51947)
@@ -0,0 +1,87 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.shaded.relief</em> creates a raster shaded relief map based on
+current resolution settings and on sun altitude, azimuth, and
+z-exaggeration values entered by the user. If no output shademap
+name is given, the new shaded relief map is named <em><input_map
+>.shade</em>. The map is assigned a grey-scale color table.
+
+<p>The parameters controlling the shading are:
+<ol>
+<li>A raster map layer to provide elevation values for the shaded
+relief map. Typically, this would be a map layer of elevation;
+however, any raster map layer can be named.
+
+<li>The <b>altitude</b> of the sun in degrees above the horizon
+(a value between 0 and 90 degrees).
+
+<li>The <b>azimuth</b> of the sun in degrees to the east of north
+(a value between 0 and 360 degrees).
+
+<li>The scaling parameter, which compensates for a different
+horizontal <b>scale</b> than vertical scale. If <b>scale</b> is a
+number, then the ewres and nsres are multiplied by that scale to
+calculate the shading. (Default=1.0 for equivalent horizontal and
+vertical scales.)
+
+<li>The <b>zmult</b> exaggeration factor that changes the apparent relief
+for the shaded relief map. This can be any positive (or negative) floating
+point value. (Default=1.0)
+
+<li>Horizontal distances are calculated in meters, using geodesic
+distances for a latitude-longitude projection. With an elevation map
+measured in feet, the <b>units</b> option can be set to automatically
+convert meters to international feet (0.3048 meters = 1 foot) or survey
+feet (1200 / 3937 meters = 1 foot). The <b>units</b> parameter overrides
+the <b>scale</b> parameter.
+
+</ol>
+
+<p><em>r.shaded.relief</em> assigns a grey-scale color table to the new
+shaded relief map.
+
+<h2>NOTES</h2>
+
+To visually improve the result of shade maps from low resolution elevation
+models, use <em>r.resamp.interp</em> with bilinear or bicubic method to
+resample the DEM at higher resolution. <em>r.shaded.relief</em> is then
+run on the resampled DEM.
+
+<h2>EXAMPLES</h2>
+
+In this example, the aspect map in the North Carolina sample
+dataset location is used to hillshade the elevation map:
+
+<div class="code"><pre>
+g.region rast=elevation -p
+r.shaded.relief input=elevation output=elevation.shaded
+</pre></div>
+
+<p>In Latitude-Longitude locations (or other non-metric locations), the
+<em>scale</em> factor has to be used:
+
+<div class="code"><pre>
+# Latitude-Longitude example
+r.shaded.relief input=srtm output=srtm.shaded scale=111120
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<p><!-- RGB version not ported to GRASS 6 (why?):
+ <em><a href="shade.clr.sh.html">shade.clr.sh</a></em><br>
+ -->
+<em>
+<a href="d.his.html">d.his</a>,
+<a href="g.region.html">g.region</a>,
+<a href="r.blend.html">r.blend</a>,
+<a href="r.colors.html">r.colors</a>,
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="r.resamp.interp.html">r.resamp.interp</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Jim Westervelt, U.S. Army Construction Engineering
+Research Laboratory
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass/trunk/raster/r.shaded.relief2/r.shaded.relief2.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Date
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list