[GRASS-SVN] r51322 - grass-addons/grass7/raster/r.convergence
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Apr 9 14:47:41 EDT 2012
Author: neteler
Date: 2012-04-09 11:47:40 -0700 (Mon, 09 Apr 2012)
New Revision: 51322
Modified:
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/slope_aspect.c
Log:
run indent on code; fix png
Modified: grass-addons/grass7/raster/r.convergence/conv.png
===================================================================
(Binary files differ)
Modified: grass-addons/grass7/raster/r.convergence/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.convergence/local_proto.h 2012-04-09 18:42:00 UTC (rev 51321)
+++ grass-addons/grass7/raster/r.convergence/local_proto.h 2012-04-09 18:47:40 UTC (rev 51322)
@@ -6,34 +6,34 @@
#include <grass/raster.h>
#ifdef MAIN
-# define GLOBAL
+#define GLOBAL
#else
-# define GLOBAL extern
+#define GLOBAL extern
#endif
/*
-PI2= PI/2
-PI4= PI/4
-M2PI=2*PI
-*/
+ PI2= PI/2
+ PI4= PI/4
+ M2PI=2*PI
+ */
#ifndef PI2
- #define PI2 (2*atan(1))
+#define PI2 (2*atan(1))
#endif
#ifndef PI4
- #define PI4 (atan(1))
+#define PI4 (atan(1))
#endif
#ifndef PI
- #define PI (4*atan(1))
+#define PI (4*atan(1))
#endif
#ifndef M2PI
- #define M2PI (8*atan(1))
+#define M2PI (8*atan(1))
#endif
#ifndef PI2PERCENT
- #define PI2PERCENT (50/atan(1))
+#define PI2PERCENT (50/atan(1))
#endif
@@ -43,43 +43,45 @@
#define RAD2DEGREE(a) ((a)*(180/PI))
-typedef char* STRING;
-typedef enum {
- m_STANDARD,
- m_INVERSE,
- m_POWER,
- m_SQUARE,
- m_GENTLE
+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 */
+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;
+typedef struct
+{
+ double cat;
+int r, g, b} FCOLORS;
-GLOBAL int gradient, f_circular, f_slope, f_method, window_size,radius;
+GLOBAL int gradient, f_circular, f_slope, f_method, window_size, radius;
GLOBAL float *aspect_matrix, *distance_matrix;
GLOBAL MAPS elevation;
-GLOBAL FCELL** slope;
-GLOBAL FCELL** aspect;
+GLOBAL FCELL **slope;
+GLOBAL FCELL **aspect;
GLOBAL int nrows, ncols;
-GLOBAL double H,V;
+GLOBAL double H, V;
GLOBAL struct Cell_head window;
-int open_map(MAPS* rast);
+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_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);
+int free_map(FCELL ** map, int n);
Modified: grass-addons/grass7/raster/r.convergence/main.c
===================================================================
--- grass-addons/grass7/raster/r.convergence/main.c 2012-04-09 18:42:00 UTC (rev 51321)
+++ grass-addons/grass7/raster/r.convergence/main.c 2012-04-09 18:47:40 UTC (rev 51322)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: r.convergence
@@ -20,166 +21,177 @@
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]);
+ struct Option *map_dem,
+ *map_slope,
+ *map_aspect, *par_window, *par_method, *par_differnce, *map_output;
- 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");
+ struct History history;
+ struct Colors colors;
+ struct Flag *flag_slope, *flag_circular;
- 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)");
-
+ 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)
+
+ 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, G_FATAL_EXIT);
- nrows = Rast_window_rows();
- ncols = Rast_window_cols();
- Rast_get_window(&window);
- radius=window_size/2;
+ 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,
+ G_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;
+ int n = 0;
float *local_aspect;
float convergence;
- out_fd = Rast_open_new(map_output->answer, FCELL_TYPE);
+ out_fd = Rast_open_new(map_output->answer, FCELL_TYPE);
out_buf = Rast_allocate_buf(FCELL_TYPE);
- strcpy(elevation.elevname,map_dem->answer);
+ 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);
+ 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);
+ 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))
+ if (row > radius && row < nrows - (radius + 1))
shift_buffers(row);
- Rast_put_row(out_fd, out_buf, FCELL_TYPE);
-
- } /* end for row */
+ Rast_put_row(out_fd, out_buf, FCELL_TYPE);
+
+ } /* end for row */
G_percent(row, nrows, 2);
-} /* end block */
+ } /* 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}};
+ {
+ 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);
-}
-
+ 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);
+ 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);
+ G_message("Done!");
+ exit(EXIT_SUCCESS);
}
Modified: grass-addons/grass7/raster/r.convergence/memory.c
===================================================================
--- grass-addons/grass7/raster/r.convergence/memory.c 2012-04-09 18:42:00 UTC (rev 51321)
+++ grass-addons/grass7/raster/r.convergence/memory.c 2012-04-09 18:47:40 UTC (rev 51322)
@@ -1,138 +1,143 @@
#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);
+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 */
+ tmp_buf = Rast_allocate_buf(rast->raster_type);
+ rast->elev = (FCELL **) G_malloc((window_size + 1) * sizeof(FCELL *));
-G_free(tmp_buf);
-return 0;
+ 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) {
+int get_cell(int col, float *buf_row, void *buf, RASTER_MAP_TYPE raster_type)
+{
- switch (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 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];
- }
-}
+ 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 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;
+ int i;
+ int col;
+ void *tmp_buf;
+ FCELL *tmp_elev_buf, *slope_tmp, *aspect_tmp;
- tmp_buf=Rast_allocate_buf(elevation.raster_type);
+ tmp_buf = Rast_allocate_buf(elevation.raster_type);
- tmp_elev_buf=elevation.elev[0];
-
- for (i = 1; i < window_size+1; ++i)
+ 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);
+ elevation.elev[window_size] = tmp_elev_buf;
+ Rast_get_row(elevation.fd, tmp_buf, row + radius + 1,
+ 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);
+ 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;
+ }
+
+ get_slope_aspect(window_size - 1);
+
+ return 0;
}
-int free_map (FCELL **map, int n) {
- int i;
- for (i=0;i<n;++i)
+int free_map(FCELL ** map, int n)
+{
+ int i;
+
+ for (i = 0; i < n; ++i)
G_free(map[i]);
- G_free(map);
- return 0;
+ G_free(map);
+ return 0;
}
-
Modified: grass-addons/grass7/raster/r.convergence/slope_aspect.c
===================================================================
--- grass-addons/grass7/raster/r.convergence/slope_aspect.c 2012-04-09 18:42:00 UTC (rev 51321)
+++ grass-addons/grass7/raster/r.convergence/slope_aspect.c 2012-04-09 18:47:40 UTC (rev 51322)
@@ -1,219 +1,256 @@
#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);
+int get_distance(int once, int row)
+{
- } else {
+ double north, south, east, west, middle;
+ double zfactor = 1;
- 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);
+ if (once) {
- }
-
- V = G_distance(east, north, east, south) / zfactor;
- H = G_distance(east, middle, west, middle) / zfactor;
+ 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);
- }
- }
+int create_distance_aspect_matrix(int row)
+{
- return 0;
-
-}
+ 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;
-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;
- }
+ aspect_matrix = G_calloc(window_size * window_size, sizeof(float));
+ distance_matrix = G_calloc(window_size * window_size, sizeof(float));
- 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);
-
+ 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;
+
+ 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;
+int get_slope_aspect(int row)
+{
- if(f_slope) {
- cur_northing=Rast_row_to_northing(row+0.5, &window);
- cur_easting=Rast_col_to_easting(col+0.5, &window);
- }
+ 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;
+ }
- 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"));
- }
+ 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);
- sum += fabs(conv)*(1/distance_modifier);
- div += 1/distance_modifier;
-
- } /* end for i, j */
+ uprow = elevation.elev[row - 1];
+ thisrow = elevation.elev[row];
+ downrow = elevation.elev[row + 1];
- sum /= div;
- sum -=PI2;
-return sum;
+ 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