[GRASS-SVN] r42493 - grass-addons/raster/r.convergence
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 7 04:09:09 EDT 2010
Author: jarekj71
Date: 2010-06-07 04:09:08 -0400 (Mon, 07 Jun 2010)
New Revision: 42493
Added:
grass-addons/raster/r.convergence/Makefile
grass-addons/raster/r.convergence/conv.png
grass-addons/raster/r.convergence/description.html
grass-addons/raster/r.convergence/local_proto.h
grass-addons/raster/r.convergence/main.c
grass-addons/raster/r.convergence/memory.c
grass-addons/raster/r.convergence/slope_aspect.c
Log:
files for r.convergence
Added: grass-addons/raster/r.convergence/Makefile
===================================================================
--- grass-addons/raster/r.convergence/Makefile (rev 0)
+++ grass-addons/raster/r.convergence/Makefile 2010-06-07 08:09:08 UTC (rev 42493)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.convergence
+
+LIBES = $(ROWIOLIB) $(GISLIB)
+DEPENDENCIES = $(ROWIODEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/raster/r.convergence/conv.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/raster/r.convergence/conv.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/raster/r.convergence/description.html
===================================================================
--- grass-addons/raster/r.convergence/description.html (rev 0)
+++ grass-addons/raster/r.convergence/description.html 2010-06-07 08:09:08 UTC (rev 42493)
@@ -0,0 +1,65 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>-s</b></DT>
+<DD>Increase convergence if slope value is high. Slope parameter radically slow down computation time, especially if window parameter is high. If slope is used addational modifier is used according to formula: sin(current)*sin(target) + cos(current)*cos(target). if slope of current and target cells are equal. The modifier is 1. If not, the modifier is applied with formula: acos(cos(convergence) * modifier) </DD>
+<DT><b>-c</b></DT>
+<DD>use circular window instead of suqare (default)</DD>
+<p>
+<DT><b>input</b></DT>
+<DD>Digital elevation model. Data can be of any type and any projection. To calculate relief convergnece, r.convergence uses real distance wchich is recalculated into cell distance, according formula: <br><code>distance_between_current_cell_and_traget_cell/distance_between_current_cell_and_nearest_neighbour_cell.</code> It is important if convergence is calculated for large areas in LatLong projecton.
+</DD>
+<p>
+<DT><b>weights</b></DT>
+<DD>Parameter describing the reduction of the impact of the cell due to its distance, where distance in cells:
+<ul>
+<li><b>standard:</b>no decay
+<li><b>inverse:</b>distance modifier is calculated as 1/x
+<li><b>power:</b>distance modifier is calculated as 1/(x*x)
+<li><b>power:</b>distance modifier is calculated as 1/(x*x)
+<li><b>gentle:</b>distance modifier is calculated as 1/((1-x)/(1+x))
+</ul>
+<p>
+<DT><b>window</b></DT>
+<DD>window size. Must be odd. For now there is no limits in window size. r.convergence uses window size instead of classical radius for compatibility with other GRASS programs.</DD>
+
+<p>
+<DT><b>output</b></DT>
+<DD>Map of convergence index. The values ranges from -100 (max divergent, real peaks and ridges) by 0 (planar areas) to 100 (max convergent, real pits and channels). Classical convergence index presented with degrees (-90 to 90)</DD>
+
+
+<h2>DESCRIPTION</h2>
+<center>
+<h3>How convergence index is calculated (3 x 3 window):<h3>
+<img src=conv.png border=1><br>
+</center>
+Convergence index is a terrain parameter which show the structure od the relief as a set of convergent areas (channels) and divergent areas (ridges). It represents the agreement of aspect direction of surrounding cells with the teoretical matrix direction. Convergence index is mean (or weighted mean if weights are used) aspect difference between real aspect and teoretical maximum divergent direction matrix representing ideal peak (see figure) minus 90 degres. So if there is maximum agreement with divergent matrix the convergence index is (0 - 90) * 10/9 = -100. If there is ideal sink (maximum convergence) the convergence index is (180 -90) * 10/9 = 100.
+Slope and aspect ere calculated internaly with the same formula as in r.slope.aspect
+Convergence index is very useful for analysis of lineamets especially represented by rigdes or chanell systems as well as valley recognition tool.
+<P>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.slope.aspect.html">r.slope.aspect</a>,
+<a href="r.param.scale.html">r.param.scale</a>,
+<a href="r.neighbour.html">r.neighbour</a>,
+</em>
+
+
+<h2>REFERENCES</h2>
+</p>
+Claps, P., Fiorentino, M., Oliveto, G., (1994), <i>Informational entropy of fractal river networks</i>,
+Journal of Hydrology, 187(1-2), 145-156 .</p>
+<p>
+Bauer J., Rohdenburg H., Bork H.-R., (1985), Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse, <i>IN: Bork, H.-R., Rohdenburg, H., Landschaftsgenese und Landschaftsoekologie, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle, Grundlagenarbeiten zu Analyse von Agrar-Oekosystemen</i>, 1-15.
+<p>
+Böhner J., Blaschke T., Montanarella, L. (eds.) (2008). SAGA Seconds Out. Hamburger Beiträge zur Physischen Geographie und Landschaftsökologie, 19: 113 s.</p>
+
+<h2>AUTHOR</h2>
+Jarek Jasiewicz
+
+<HR>
+<P><a href="index.html">Main index</a> - <a href="raster.html">raster index</a> - <a href="full_index.html">Full index</a></P>
+<P>© 2003-2009 <a href="http://grass.osgeo.org">GRASS Development Team</a></p>
+</body>
+</html>
Added: grass-addons/raster/r.convergence/local_proto.h
===================================================================
--- grass-addons/raster/r.convergence/local_proto.h (rev 0)
+++ grass-addons/raster/r.convergence/local_proto.h 2010-06-07 08:09:08 UTC (rev 42493)
@@ -0,0 +1,80 @@
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+/*
+PI2= PI/2
+PI4= PI/4
+M2PI=2*PI
+*/
+#ifndef PI2
+ #define PI2 (2*atan(1))
+#endif
+
+#ifndef PI4
+ #define PI4 (atan(1))
+#endif
+
+#ifndef PI
+ #define PI (4*atan(1))
+#endif
+
+#ifndef M2PI
+ #define M2PI (8*atan(1))
+#endif
+
+#ifndef PI2PERCENT
+ #define PI2PERCENT (50/atan(1))
+#endif
+
+
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define DEGREE2RAD(a) ((a)/(180/PI))
+#define RAD2DEGREE(a) ((a)*(180/PI))
+
+
+typedef char* STRING;
+typedef enum {
+ m_STANDARD,
+ m_INVERSE,
+ m_POWER,
+ m_SQUARE,
+ m_GENTLE
+} methods;
+
+typedef struct {
+ char elevname[150];
+ RASTER_MAP_TYPE raster_type;
+ FCELL **elev;
+ int fd; /* file descriptor */
+} MAPS;
+
+int gradient, f_circular, f_slope, f_method, window_size,radius;
+float *aspect_matrix, *distance_matrix;
+MAPS elevation;
+FCELL** slope;
+FCELL** aspect;
+
+int nrows, ncols;
+double H,V;
+struct Cell_head window;
+
+int open_map(MAPS* rast);
+int create_maps(void);
+int shift_buffers(int row);
+int get_cell (int col, float* buf_row, void* buf, RASTER_MAP_TYPE raster_type);
+int get_slope_aspect(int row);
+int get_distance(int once, int row);
+int create_distance_aspect_matrix(int row);
+float calculate_convergence(int row, int cur_row, int col);
+int free_map(FCELL **map, int n);
Added: grass-addons/raster/r.convergence/main.c
===================================================================
--- grass-addons/raster/r.convergence/main.c (rev 0)
+++ grass-addons/raster/r.convergence/main.c 2010-06-07 08:09:08 UTC (rev 42493)
@@ -0,0 +1,182 @@
+/****************************************************************************
+ *
+ * MODULE: r.convergence
+ * AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl
+ * Original convergence index in SAGA GIS sofware: Olaf Conrad
+ *
+ * PURPOSE: Calculate convergence index (parameter defining the local convergence of the relief)
+*
+*
+* COPYRIGHT: (C) 2002,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.
+*
+ *****************************************************************************/
+
+#define MAIN
+#include "local_proto.h"
+
+int main(int argc, char **argv)
+{
+ struct Option *map_dem,
+ *map_slope,
+ *map_aspect,
+ *par_window,
+ *par_method,
+ *par_differnce,
+ *map_output;
+
+ struct History history;
+ struct Colors colors;
+ struct Flag *flag_slope, *flag_circular;
+
+ int outfd;
+ float tmp;
+ FCELL *out_buf;
+
+ int i,j, n;
+ G_gisinit(argv[0]);
+
+ map_dem = G_define_standard_option(G_OPT_R_INPUT);
+ map_dem->description = _("Digital elevation model map");
+
+ map_output = G_define_standard_option(G_OPT_R_OUTPUT);
+ map_output->description = _("Output convergence index map");
+
+ par_window = G_define_option();
+ par_window->key = "window";
+ par_window->type = TYPE_INTEGER;
+ par_window->answer = "3";
+ par_window->required = YES;
+ par_window->description = _("Window size");
+
+ par_method = G_define_option();
+ par_method->key = "weights";
+ par_method->type = TYPE_STRING;
+ par_method->options = "standard,inverse,power,square,gentle";
+ par_method->answer = "standard";
+ par_method->required = YES;
+ par_method->description = _("Method for reducing the impact of the cell due to distance");
+
+ flag_circular = G_define_flag();
+ flag_circular->key = 'c';
+ flag_circular->description = _("Use circular window (default: square)");
+
+ flag_slope = G_define_flag();
+ flag_slope->key = 's';
+ flag_slope->description = _("Add slope convergence (radically slow down calculation time)");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ window_size=atoi(par_window->answer);
+ if(window_size<3 || window_size%2==0)
+ G_fatal_error(_("Window size must be odd and at least 3"));
+
+ if (!strcmp(par_method->answer, "standard"))
+ f_method=m_STANDARD;
+ else if (!strcmp(par_method->answer, "inverse"))
+ f_method=m_INVERSE;
+ else if (!strcmp(par_method->answer, "power"))
+ f_method=m_POWER;
+ else if (!strcmp(par_method->answer, "square"))
+ f_method=m_SQUARE;
+ else if (!strcmp(par_method->answer, "gentle"))
+ f_method=m_GENTLE;
+
+ f_circular=(flag_circular->answer !=0);
+ f_slope=(flag_slope->answer !=0);
+
+ /* align window */
+ G_check_input_output_name(map_dem->answer, map_output->answer, GR_FATAL_EXIT);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ radius=window_size/2;
+
+
+{
+ int row, col;
+ int cur_row;
+ int i_row, i_col, j_row, j_col;
+ float contr_aspect;
+ float dir;
+ int n=0;
+ float *local_aspect;
+ float convergence;
+
+ strcpy(elevation.elevname,map_dem->answer);
+ open_map(&elevation);
+ create_maps();
+
+ /* create aspect matrix and distance_matrix */
+ create_distance_aspect_matrix(0);
+
+ if ((outfd = G_open_raster_new(map_output->answer, FCELL_TYPE)) < 0)
+ G_fatal_error(_("Unable to create raster map <%s>"), map_output->answer);
+
+ out_buf = G_allocate_raster_buf(FCELL_TYPE);
+
+
+ /* cur_row: current to row in the buffer */
+ /* open_map and create_maps create data for first pass */
+
+
+ for(row=0;row<nrows;++row) {
+ G_percent(row, nrows, 2);
+ cur_row = (row < radius ) ? row :
+ ((row >= nrows-radius) ? window_size - (nrows-row) : radius);
+
+ if (row<1 || row >nrows-1) {
+ G_set_f_null_value(out_buf,ncols);
+ continue;
+ }
+
+ for (col=0;col<ncols;++col) {
+ if (col<1 || col >ncols-1) {
+ G_set_f_null_value(&out_buf[col],1);
+ continue;
+ }
+ out_buf[col]=PI2PERCENT*calculate_convergence(row,cur_row,col);
+ }
+
+ if(row>radius && row<nrows-(radius+1)) {
+ shift_buffers(row);
+ }
+
+ if (G_put_raster_row(outfd, out_buf, FCELL_TYPE) < 0)
+ G_fatal_error(_("Failed writing raster map <%s>"), map_output->answer);
+ } /* end for row */
+ G_percent(row, nrows, 2);
+} /* end block */
+
+
+ G_init_colors(&colors);
+ G_add_color_rule(-100, 56, 0, 0, -70, 128, 0, 0, &colors);
+ G_add_color_rule(-70, 128, 0, 0, -50, 255, 0, 0, &colors);
+ G_add_color_rule(-50, 255, 0, 0, 0, 255, 255, 255, &colors);
+ G_add_color_rule(0, 255, 255, 255, 50, 0, 0, 255, &colors);
+ G_add_color_rule(50, 0, 0, 255, 70, 0, 0, 128, &colors);
+ G_add_color_rule(70, 0, 0, 128, 100, 0, 0, 56, &colors);
+ G_write_colors(map_output->answer, G_mapset(), &colors);
+
+
+ free_map(slope, window_size);
+ free_map(aspect, window_size);
+ free_map(elevation.elev, window_size+1);
+ G_free(distance_matrix);
+ G_free(aspect_matrix);
+
+ G_free(out_buf);
+ G_close_cell(outfd);
+
+ G_short_history(map_output->answer, "raster", &history);
+ G_command_history(&history);
+ G_write_history(map_output->answer, &history);
+
+
+G_message("Done!");
+exit(EXIT_SUCCESS);
+
+}
Added: grass-addons/raster/r.convergence/memory.c
===================================================================
--- grass-addons/raster/r.convergence/memory.c (rev 0)
+++ grass-addons/raster/r.convergence/memory.c 2010-06-07 08:09:08 UTC (rev 42493)
@@ -0,0 +1,150 @@
+#include "local_proto.h"
+
+int open_map(MAPS* rast) {
+
+ int row, col;
+ int fd;
+ char* mapset;
+ struct Cell_head cellhd;
+ int bufsize;
+ void* tmp_buf;
+
+ mapset = G_find_cell2(rast->elevname, "");
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), rast->elevname);
+
+ if ((rast->fd = G_open_cell_old(rast->elevname, mapset)) < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), rast->elevname);
+
+ if (G_get_cellhd(rast->elevname, mapset, &cellhd) < 0)
+ G_fatal_error(_("Unable to read file header of <%s>"), rast->elevname);
+
+ rast->raster_type = G_get_raster_map_type(rast->fd);
+
+ tmp_buf=G_allocate_raster_buf(rast->raster_type);
+ rast->elev = (FCELL**) G_malloc((window_size+1) * sizeof(FCELL*));
+
+ for (row = 0; row < window_size+1; ++row) {
+ rast->elev[row] = G_allocate_raster_buf(FCELL_TYPE);
+
+ if (G_get_raster_row(rast->fd, tmp_buf,row, rast->raster_type)<0) {
+ G_close_cell(rast->fd);
+ G_fatal_error(_("Cannot to read <%s> at row <%d>"), rast->elevname,row);
+ }
+ for (col=0;col<ncols;++col)
+ get_cell(col, rast->elev[row], tmp_buf, rast->raster_type);
+ } /* end elev */
+
+G_free(tmp_buf);
+return 0;
+}
+
+int get_cell(int col, float* buf_row, void* buf, RASTER_MAP_TYPE raster_type) {
+
+ switch (raster_type) {
+
+ case CELL_TYPE:
+ if (G_is_null_value(&((CELL *) buf)[col],CELL_TYPE))
+ G_set_f_null_value(&buf_row[col],1);
+ else
+ buf_row[col] = (FCELL) ((CELL *) buf)[col];
+ break;
+
+ case FCELL_TYPE:
+ if (G_is_null_value(&((FCELL *) buf)[col],FCELL_TYPE))
+ G_set_f_null_value(&buf_row[col],1);
+ else
+ buf_row[col] = (FCELL) ((FCELL *) buf)[col];
+ break;
+
+ case DCELL_TYPE:
+ if (G_is_null_value(&((DCELL *) buf)[col],DCELL_TYPE))
+ G_set_f_null_value(&buf_row[col],1);
+ else
+ buf_row[col] = (FCELL) ((DCELL *) buf)[col];
+ break;
+ }
+
+return 0;
+}
+
+
+int create_maps(void)
+{
+ int row,col;
+
+ G_begin_distance_calculations();
+ if(G_projection() != PROJECTION_LL)
+ get_distance(1,0);
+
+ slope = (FCELL**) G_malloc(window_size * sizeof(FCELL*));
+ aspect = (FCELL**) G_malloc(window_size * sizeof(FCELL*));
+ for (row = 0; row < window_size; ++row) {
+
+ if(G_projection() == PROJECTION_LL) {
+ get_distance(0,row);
+ create_distance_aspect_matrix(row);
+ }
+
+ slope[row] = G_allocate_raster_buf(FCELL_TYPE);
+ aspect[row] = G_allocate_raster_buf(FCELL_TYPE);
+ get_slope_aspect(row);
+ }
+}
+
+int shift_buffers(int row)
+{
+ int i;
+ int col;
+ void* tmp_buf;
+ FCELL* tmp_elev_buf, *slope_tmp, *aspect_tmp;
+
+ tmp_buf=G_allocate_raster_buf(elevation.raster_type);
+
+ tmp_elev_buf=elevation.elev[0];
+
+ for (i = 1; i < window_size+1; ++i)
+ elevation.elev[i - 1] = elevation.elev[i];
+
+ elevation.elev[window_size]=tmp_elev_buf;
+
+ if (G_get_raster_row(elevation.fd, tmp_buf,row+radius+1, elevation.raster_type)<0) {
+ G_close_cell(elevation.fd);
+ G_fatal_error(_("Cannot read <%s> at row <%d>"), elevation.elevname,row);
+ }
+
+ for (col=0;col<ncols;++col)
+ get_cell(col, elevation.elev[window_size], tmp_buf, elevation.raster_type);
+
+ G_free(tmp_buf);
+
+ slope_tmp = slope[0];
+ aspect_tmp = aspect[0];
+
+ for (i = 1; i < window_size; ++i) {
+ slope[i - 1] = slope[i];
+ aspect[i - 1] = aspect[i];
+ }
+
+ slope[window_size-1] = slope_tmp;
+ aspect[window_size-1] = aspect_tmp;
+
+ if(G_projection() == PROJECTION_LL) {
+ get_distance(0,row);
+ create_distance_aspect_matrix(row);
+ }
+
+ get_slope_aspect(window_size-1);
+
+ return 0;
+}
+
+int free_map (FCELL **map, int n) {
+ int i;
+ for (i=0;i<n;++i)
+ G_free(map[i]);
+ G_free(map);
+ return 0;
+}
+
Added: grass-addons/raster/r.convergence/slope_aspect.c
===================================================================
--- grass-addons/raster/r.convergence/slope_aspect.c (rev 0)
+++ grass-addons/raster/r.convergence/slope_aspect.c 2010-06-07 08:09:08 UTC (rev 42493)
@@ -0,0 +1,230 @@
+#include "local_proto.h"
+int get_distance(int once, int row) {
+
+ double north,south,east,west,middle;
+ double zfactor=1;
+
+ if(once) {
+
+ G_get_window(&window);
+ north = G_row_to_northing(0.5, &window);
+ middle = G_row_to_northing(1.5, &window);
+ south = G_row_to_northing(2.5, &window);
+ east = G_col_to_easting(2.5, &window);
+ west = G_col_to_easting(0.5, &window);
+
+ } else {
+
+ north = G_row_to_northing(row+0.5, &window);
+ middle = G_row_to_northing(row+1.5, &window);
+ south = G_row_to_northing(row+2.5, &window);
+ east = G_col_to_easting(2.5, &window);
+ west = G_col_to_easting(0.5, &window);
+
+ }
+
+ V = G_distance(east, north, east, south) / zfactor;
+ H = G_distance(east, middle, west, middle) / zfactor;
+ return 0;
+
+}
+
+int create_distance_aspect_matrix(int row) {
+
+ int i;
+ int mx,x,my,y;
+ int dx,dy;
+ float distance;
+ int set_to_zero=0;
+ float cur_northing, cur_easting, target_northing, target_easting;
+ float ns_dist, ew_dist;
+ float min_cell_size;
+
+ aspect_matrix=G_calloc(window_size*window_size,sizeof(float));
+ distance_matrix=G_calloc(window_size*window_size,sizeof(float));
+
+ cur_northing=G_row_to_northing(row+0.5, &window);
+ cur_easting=G_col_to_easting(0.5, &window);
+ target_northing=G_row_to_northing(row+1.5, &window);
+ target_easting=G_col_to_easting(1.5, &window);
+
+ ns_dist=G_distance(cur_easting, cur_northing, cur_easting,target_northing);
+ ew_dist=G_distance(cur_easting, cur_northing, target_easting,cur_northing);
+
+ min_cell_size=MIN(ns_dist,ew_dist);
+
+ for(my=0, y=-radius;my<window_size;++my, ++y)
+ for(mx=0, x=radius;mx<window_size;++mx, --x) {
+
+ cur_northing=G_row_to_northing(row+0.5, &window);
+ cur_easting=G_col_to_easting(0.5, &window);
+ target_northing=G_row_to_northing(row+y+0.5, &window);
+ target_easting=G_col_to_easting(x+0.5, &window);
+
+ distance=G_distance(cur_easting, cur_northing, target_easting,target_northing);
+ distance /= min_cell_size;
+
+ set_to_zero=0;
+ if(distance<1)
+ set_to_zero=1;
+ if(f_circular && distance > radius)
+ set_to_zero=1;
+
+
+ if(set_to_zero)
+ {
+ distance_matrix[my*window_size+mx]=0;
+ aspect_matrix[my*window_size+mx]=-1;
+ }
+ else
+ {
+ distance_matrix[my*window_size+mx]=distance;
+ ns_dist=G_distance(cur_easting, cur_northing, cur_easting,target_northing);
+ ns_dist=(cur_northing>target_northing) ? ns_dist : -ns_dist;
+ ew_dist=G_distance(cur_easting, cur_northing, target_easting,cur_northing);
+ ew_dist=(cur_easting>target_easting) ? -ew_dist : ew_dist;
+
+ aspect_matrix[my*window_size+mx]=
+ (y != 0.) ? PI + atan2(ew_dist, ns_dist) :
+ (x>0. ? PI2+PI : PI2);
+ }
+ }
+
+ return 0;
+
+}
+
+int get_slope_aspect(int row) {
+
+ int col, i;
+ FCELL c[10];
+ FCELL dx, dy;
+ FCELL *uprow,*thisrow,*downrow;
+ int d_row[]={-1,-1 ,-1, 0, 0, 0, 1, 1, 1};
+ int d_col[]={-1, 0 ,1, -1,0 ,1, -1, 0 ,1};
+
+ if(row<1 || row > nrows-2) {
+ G_set_f_null_value(slope[row],ncols);
+ G_set_f_null_value(aspect[row],ncols);
+ return 1;
+ }
+
+ G_set_f_null_value(&slope[row][0],1);
+ G_set_f_null_value(&aspect[row][0],1);
+ G_set_f_null_value(&slope[row][ncols-1],1);
+ G_set_f_null_value(&aspect[row][ncols-1],1);
+
+ uprow=elevation.elev[row-1];
+ thisrow=elevation.elev[row];
+ downrow=elevation.elev[row+1];
+
+ /*
+ c1 = uprow
+ c2 = uprow+1
+ c3 = uprow+2
+ c4 = thisrow
+ c5 = thisrow+1
+ c6 = thisrow+2
+ c7 = downrow
+ c8 = downrow+1
+ c9 = downrow+2
+ */
+
+
+ for (col=1;col<ncols-1;++col) {
+
+ for (i=0;i<9;++i)
+ if(G_is_f_null_value(&elevation.elev[row+d_row[i]][col+d_col[i]])) {
+ G_set_f_null_value(&slope[row][col],1);
+ G_set_f_null_value(&aspect[row][col],1);
+ continue;
+ }
+
+ dx= ((elevation.elev[row-1][col-1] + 2 * elevation.elev[row][col-1] + elevation.elev[row+1][col-1]) -
+ (elevation.elev[row-1][col+1] + 2 * elevation.elev[row][col+1] + elevation.elev[row+1][col+1])) /(H*4);
+
+ dy= ((elevation.elev[row+1][col-1] + 2 * elevation.elev[row+1][col] + elevation.elev[row+1][col+1]) -
+ (elevation.elev[row-1][col-1] + 2 * elevation.elev[row-1][col] + elevation.elev[row-1][col+1])) /(V*4);
+
+ slope[row][col] = atan(sqrt(dx * dx + dy * dy));
+ if(dy==0.)
+ aspect[row][col] = dx<0 ? PI+PI2 : PI2;
+ else
+ aspect[row][col] = atan2(dx, dy)<0 ? atan2(dx, dy) + M2PI : atan2(dx, dy);
+
+ }
+ return 0;
+}
+
+float calculate_convergence(int row, int cur_row, int col) {
+ int i, j, k=0;
+ float i_distance;
+ float conv, sum=0, div=0;
+ float x;
+ float cur_slope;
+ float slope_modifier;
+ float distance_modifier;
+ float cur_northing, cur_easting, target_northing, target_easting, cur_distance;
+
+ for(i=-radius;i < radius+1;++i)
+ for(j=-radius;j < radius+1;++j, ++k) {
+
+ if(cur_row+i<0 || cur_row+i>window_size-1 || col+j<1 || col +j > ncols-2)
+ continue;
+ x=distance_matrix[k];
+
+ if(x<1)
+ continue;
+
+ if(f_slope) {
+
+ cur_northing=G_row_to_northing(row+0.5, &window);
+ cur_easting=G_col_to_easting(col+0.5, &window);
+ target_northing=G_row_to_northing(row+i+0.5, &window);
+ target_easting=G_col_to_easting(col+j+0.5, &window);
+
+ cur_distance=G_distance(cur_easting, cur_northing, target_easting,target_northing);
+ cur_slope=(atan((elevation.elev[cur_row+i][col+j]-elevation.elev[cur_row][col])/cur_distance));
+ slope_modifier = sin(slope[cur_row][col])*sin(cur_slope)+
+ cos(slope[cur_row][col])*cos(cur_slope);
+ }
+
+ conv = aspect[cur_row+i][col+j]-aspect_matrix[k];
+
+ if(f_slope)
+ conv = acos(slope_modifier*cos(conv));
+
+ if (conv < -PI)
+ conv += M2PI;
+ else if (conv > PI)
+ conv -= M2PI;
+
+ switch (f_method) {
+ case 0:
+ distance_modifier=1;
+ break;
+ case 1: /* inverse */
+ distance_modifier=x;
+ break;
+ case 2: /* power */
+ distance_modifier=x*x;
+ break;
+ case 3: /* square */
+ distance_modifier=sqrt(x);
+ break;
+ case 4: /* gentle */
+ distance_modifier=1+((1-x)*(1+x));
+ break;
+ default:
+ G_fatal_error(_("Decay: wrong option"));
+ }
+
+ sum += fabs(conv)*(1/distance_modifier);
+ div += 1/distance_modifier;
+
+ } /* end for i, j */
+
+ sum /= div;
+ sum -=PI2;
+return sum;
+}
More information about the grass-commit
mailing list