/**************************************************************************** * * MODULE: r.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 * 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 * 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 * 11/99 updated $ewres to ewres() and $nsres to nsres() * updated number to FP in r.mapcalc statement Markus Neteler */ #include #include #include #include #include #include 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; int overwrite=0; 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(_("elevation")); G_add_keyword(_("relief")); G_add_keyword(_("terrain")); G_add_keyword(_("hillshade")); 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->label = _("Name for output shaded relief map"); 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.altitude->guisection = _("Sun position"); 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.azimuth->guisection = _("Sun position"); parm.zmult = G_define_option(); parm.zmult->key = "zscale"; 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; /* Check due to default output map = input.shade map */ overwrite = G_check_overwrite(argc, argv); 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 (G_find_raster2(sr_name, G_mapset())) { if (overwrite) G_warning(_("Raster map <%s> already exists and will be overwritten"), sr_name); else G_fatal_error(_("Raster map <%s> already exists"), sr_name); } 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.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); }