[GRASS-SVN] r45420 - in grass-addons/grass7/raster: . r.convergence
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 17 14:40:58 EST 2011
Author: jarekj71
Date: 2011-02-17 11:40:58 -0800 (Thu, 17 Feb 2011)
New Revision: 45420
Added:
grass-addons/grass7/raster/r.convergence/
grass-addons/grass7/raster/r.convergence/Makefile
grass-addons/grass7/raster/r.convergence/conv.png
grass-addons/grass7/raster/r.convergence/local_proto.h
grass-addons/grass7/raster/r.convergence/main.c
grass-addons/grass7/raster/r.convergence/memory.c
grass-addons/grass7/raster/r.convergence/r.convergence.html
grass-addons/grass7/raster/r.convergence/r.convergence.tmp.html
grass-addons/grass7/raster/r.convergence/slope_aspect.c
Log:
Initial import
Added: grass-addons/grass7/raster/r.convergence/Makefile
===================================================================
--- grass-addons/grass7/raster/r.convergence/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.convergence/Makefile 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.convergence
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.convergence/conv.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.convergence/conv.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.convergence/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.convergence/local_proto.h (rev 0)
+++ grass-addons/grass7/raster/r.convergence/local_proto.h 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,85 @@
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include <grass/gis.h>
+#include <grass/raster.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;
+
+typedef struct {
+ double cat;
+ int r,g,b
+} FCOLORS;
+
+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/grass7/raster/r.convergence/main.c
===================================================================
--- grass-addons/grass7/raster/r.convergence/main.c (rev 0)
+++ grass-addons/grass7/raster/r.convergence/main.c 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,185 @@
+/****************************************************************************
+ *
+ * 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 out_fd;
+ 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 = Rast_window_rows();
+ ncols = Rast_window_cols();
+ Rast_get_window(&window);
+ 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;
+
+ out_fd = Rast_open_new(map_output->answer, FCELL_TYPE);
+ out_buf = Rast_allocate_buf(FCELL_TYPE);
+
+ strcpy(elevation.elevname,map_dem->answer);
+ open_map(&elevation);
+ create_maps();
+
+ /* create aspect matrix and distance_matrix */
+ create_distance_aspect_matrix(0);
+
+ /* 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<radius || row >nrows-radius) {
+ Rast_set_f_null_value(out_buf,ncols);
+ Rast_put_row(out_fd, out_buf, FCELL_TYPE);
+ continue;
+ }
+
+ for (col=0;col<ncols;++col) {
+ if (col<radius || col >ncols-radius) {
+ Rast_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);
+
+ Rast_put_row(out_fd, out_buf, FCELL_TYPE);
+
+ } /* end for row */
+ G_percent(row, nrows, 2);
+} /* end block */
+
+{
+ FCOLORS fcolr[7]={
+ {-100, 56, 0, 0},
+ {-70, 128, 0, 0},
+ {-50, 255, 0, 0},
+ {0, 255, 255, 255},
+ {50, 0, 0, 255},
+ {70, 0, 0, 128},
+ {100, 0, 0, 56}};
+
+ Rast_init_colors(&colors);
+ for(i=0;i<6;++i)
+ Rast_add_d_color_rule(&fcolr[i].cat, fcolr[i].r, fcolr[i].g, fcolr[i].b,
+ &fcolr[i+1].cat, fcolr[i+1].r, fcolr[i+1].g, fcolr[i+1].b, &colors);
+ Rast_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);
+ Rast_close(out_fd);
+
+ Rast_short_history(map_output->answer, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(map_output->answer, &history);
+
+
+G_message("Done!");
+exit(EXIT_SUCCESS);
+
+}
Added: grass-addons/grass7/raster/r.convergence/memory.c
===================================================================
--- grass-addons/grass7/raster/r.convergence/memory.c (rev 0)
+++ grass-addons/grass7/raster/r.convergence/memory.c 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,138 @@
+#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 = (char*)G_find_raster2(rast->elevname, "");
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), rast->elevname);
+
+ rast->fd = Rast_open_old(rast->elevname, mapset);
+ Rast_get_cellhd(rast->elevname, mapset, &cellhd);
+ rast->raster_type = Rast_map_type(rast->elevname, mapset);
+
+ if (window.ew_res < cellhd.ew_res || window.ns_res < cellhd.ns_res)
+ G_fatal_error(_("Region resolution shoudn't be lesser than map %s resolution. Run g.region rast=%s to set proper resolution"),
+ rast->elevname, rast->elevname);
+
+ tmp_buf=Rast_allocate_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] = Rast_allocate_buf(FCELL_TYPE);
+ Rast_get_row(rast->fd, tmp_buf,row, rast->raster_type);
+ 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 (Rast_is_null_value(&((CELL *) buf)[col],CELL_TYPE))
+ Rast_set_f_null_value(&buf_row[col],1);
+ else
+ buf_row[col] = (FCELL) ((CELL *) buf)[col];
+ break;
+
+ case FCELL_TYPE:
+ if (Rast_is_null_value(&((FCELL *) buf)[col],FCELL_TYPE))
+ Rast_set_f_null_value(&buf_row[col],1);
+ else
+ buf_row[col] = (FCELL) ((FCELL *) buf)[col];
+ break;
+
+ case DCELL_TYPE:
+ if (Rast_is_null_value(&((DCELL *) buf)[col],DCELL_TYPE))
+ Rast_set_f_null_value(&buf_row[col],1);
+ else
+ buf_row[col] = (FCELL) ((DCELL *) buf)[col];
+ }
+}
+
+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] = Rast_allocate_buf(FCELL_TYPE);
+ aspect[row] = Rast_allocate_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=Rast_allocate_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;
+ Rast_get_row(elevation.fd, tmp_buf,row+radius+1, elevation.raster_type);
+
+ 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/grass7/raster/r.convergence/r.convergence.html
===================================================================
--- grass-addons/grass7/raster/r.convergence/r.convergence.html (rev 0)
+++ grass-addons/grass7/raster/r.convergence/r.convergence.html 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,59 @@
+<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>
+
+<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>
+
+<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>
+
+<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>
+
+
+<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.
+
+
+<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
Added: grass-addons/grass7/raster/r.convergence/r.convergence.tmp.html
===================================================================
--- grass-addons/grass7/raster/r.convergence/r.convergence.tmp.html (rev 0)
+++ grass-addons/grass7/raster/r.convergence/r.convergence.tmp.html 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,52 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS manual: r.convergence</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>r.convergence</b></em> <h2>KEYWORDS</h2>
+<h2>SYNOPSIS</h2>
+<b>r.convergence</b><br>
+<b>r.convergence help</b><br>
+<b>r.convergence</b> [-<b>cs</b>] <b>input</b>=<em>name</em> <b>output</b>=<em>name</em> <b>window</b>=<em>integer</em> <b>weights</b>=<em>string</em> [--<b>overwrite</b>] [--<b>verbose</b>] [--<b>quiet</b>]
+
+<h3>Flags:</h3>
+<DL>
+<DT><b>-c</b></DT>
+<DD>Use circular window (default: square)</DD>
+
+<DT><b>-s</b></DT>
+<DD>Add slope convergence (radically slow down calculation time)</DD>
+
+<DT><b>--overwrite</b></DT>
+<DD>Allow output files to overwrite existing files</DD>
+<DT><b>--verbose</b></DT>
+<DD>Verbose module output</DD>
+<DT><b>--quiet</b></DT>
+<DD>Quiet module output</DD>
+</DL>
+
+<h3>Parameters:</h3>
+<DL>
+<DT><b>input</b>=<em>name</em></DT>
+<DD>Digital elevation model map</DD>
+
+<DT><b>output</b>=<em>name</em></DT>
+<DD>Output convergence index map</DD>
+
+<DT><b>window</b>=<em>integer</em></DT>
+<DD>Window size</DD>
+<DD>Default: <em>3</em></DD>
+
+<DT><b>weights</b>=<em>string</em></DT>
+<DD>Method for reducing the impact of the cell due to distance</DD>
+<DD>Options: <em>standard,inverse,power,square,gentle</em></DD>
+<DD>Default: <em>standard</em></DD>
+
+</DL>
Added: grass-addons/grass7/raster/r.convergence/slope_aspect.c
===================================================================
--- grass-addons/grass7/raster/r.convergence/slope_aspect.c (rev 0)
+++ grass-addons/grass7/raster/r.convergence/slope_aspect.c 2011-02-17 19:40:58 UTC (rev 45420)
@@ -0,0 +1,219 @@
+#include "local_proto.h"
+int get_distance(int once, int row) {
+
+ double north,south,east,west,middle;
+ double zfactor=1;
+
+ if(once) {
+
+ north = Rast_row_to_northing(0.5, &window);
+ middle = 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);
+
+ } else {
+
+ north = Rast_row_to_northing(row+0.5, &window);
+ middle = Rast_row_to_northing(row+1.5, &window);
+ south = Rast_row_to_northing(row+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) / 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=Rast_row_to_northing(row+0.5, &window);
+ cur_easting=Rast_col_to_easting(0.5, &window);
+ target_northing=Rast_row_to_northing(row+1.5, &window);
+ target_easting=Rast_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=Rast_row_to_northing(row+0.5, &window);
+ cur_easting=Rast_col_to_easting(0.5, &window);
+ target_northing=Rast_row_to_northing(row+y+0.5, &window);
+ target_easting=Rast_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) {
+ Rast_set_f_null_value(slope[row],ncols);
+ Rast_set_f_null_value(aspect[row],ncols);
+ return 1;
+ }
+
+ Rast_set_f_null_value(&slope[row][0],1);
+ Rast_set_f_null_value(&aspect[row][0],1);
+ Rast_set_f_null_value(&slope[row][ncols-1],1);
+ Rast_set_f_null_value(&aspect[row][ncols-1],1);
+
+ uprow=elevation.elev[row-1];
+ thisrow=elevation.elev[row];
+ downrow=elevation.elev[row+1];
+
+ for (col=1;col<ncols-1;++col) {
+
+ for (i=0;i<9;++i)
+ if(Rast_is_f_null_value(&elevation.elev[row+d_row[i]][col+d_col[i]])) {
+ Rast_set_f_null_value(&slope[row][col],1);
+ Rast_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;
+
+ if(f_slope) {
+ cur_northing=Rast_row_to_northing(row+0.5, &window);
+ cur_easting=Rast_col_to_easting(col+0.5, &window);
+ }
+
+
+ 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) {
+ target_northing=Rast_row_to_northing(row+i+0.5, &window);
+ target_easting=Rast_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