[GRASS-SVN] r43677 - in grass/trunk/raster: . r.regression.line
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Sep 24 11:24:56 EDT 2010
Author: mmetz
Date: 2010-09-24 15:24:55 +0000 (Fri, 24 Sep 2010)
New Revision: 43677
Added:
grass/trunk/raster/r.regression.line/
grass/trunk/raster/r.regression.line/Makefile
grass/trunk/raster/r.regression.line/main.c
grass/trunk/raster/r.regression.line/r.regression.line.html
Modified:
grass/trunk/raster/Makefile
Log:
add C version of r.regression.line
Modified: grass/trunk/raster/Makefile
===================================================================
--- grass/trunk/raster/Makefile 2010-09-24 15:20:38 UTC (rev 43676)
+++ grass/trunk/raster/Makefile 2010-09-24 15:24:55 UTC (rev 43677)
@@ -73,6 +73,7 @@
r.reclass \
r.recode \
r.region \
+ r.regression.line \
r.report \
r.resamp.bspline \
r.resamp.filter \
Added: grass/trunk/raster/r.regression.line/Makefile
===================================================================
--- grass/trunk/raster/r.regression.line/Makefile (rev 0)
+++ grass/trunk/raster/r.regression.line/Makefile 2010-09-24 15:24:55 UTC (rev 43677)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.regression.line
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass/trunk/raster/r.regression.line/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass/trunk/raster/r.regression.line/main.c
===================================================================
--- grass/trunk/raster/r.regression.line/main.c (rev 0)
+++ grass/trunk/raster/r.regression.line/main.c 2010-09-24 15:24:55 UTC (rev 43677)
@@ -0,0 +1,150 @@
+/*
+ * r.regression
+ *
+ * Calculates linear regression from two raster maps: y = a + b*x
+ * Copyright (C) 2010 by the GRASS Development Team
+ * Author(s): Markus Metz
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ * This program is a replacement for the r.regression.line script
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+
+int main(int argc, char *argv[])
+{
+ unsigned int r, c, rows, cols; /* totals */
+ int map1_fd, map2_fd;
+ double sumX, sumY, sumsqX, sumsqY, sumXY;
+ double meanX, meanY, varX, varY, sdX, sdY;
+ double A, B, R, F;
+ long count = 0;
+ DCELL *map1_buf, *map2_buf, map1_val, map2_val;
+ char *name;
+ struct Option *input_map1, *input_map2, *output_opt;
+ struct Flag *shell_style;
+ struct Cell_head region;
+ struct GModule *module;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("statistics"));
+ module->description =
+ _("Calculates linear regression from two raster maps: y = a + b*x.");
+
+ /* Define the different options */
+ input_map1 = G_define_standard_option(G_OPT_R_MAP);
+ input_map1->key = "map1";
+ input_map1->description = (_("Map for x coefficient"));
+
+ input_map2 = G_define_standard_option(G_OPT_R_MAP);
+ input_map2->key = "map2";
+ input_map2->description = (_("Map for y coefficient"));
+
+ output_opt = G_define_standard_option(G_OPT_F_OUTPUT);
+ output_opt->key = "output";
+ output_opt->required = NO;
+ output_opt->description =
+ (_("ASCII file for storing regression coefficients (output to screen if file not specified)."));
+
+ shell_style = G_define_flag();
+ shell_style->key = 'g';
+ shell_style->description = _("Print in shell script style");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ name = output_opt->answer;
+ if (name != NULL && strcmp(name, "-") != 0) {
+ if (NULL == freopen(name, "w", stdout)) {
+ G_fatal_error(_("Unable to open file <%s> for writing"), name);
+ }
+ }
+
+ G_get_window(®ion);
+ rows = region.rows;
+ cols = region.cols;
+
+ /* open maps */
+ map1_fd = Rast_open_old(input_map1->answer, G_mapset());
+ map2_fd = Rast_open_old(input_map2->answer, G_mapset());
+
+ map1_buf = Rast_allocate_d_buf();
+ map2_buf = Rast_allocate_d_buf();
+
+ sumX = sumY = sumsqX = sumsqY = sumXY = 0.0;
+ meanX = meanY = varX = varY = sdX = sdY = 0.0;
+ for (r = 0; r < rows; r++) {
+ G_percent(r, rows, 2);
+ Rast_get_d_row(map1_fd, map1_buf, r);
+ Rast_get_d_row(map2_fd, map2_buf, r);
+ for (c = 0; c < cols; c++) {
+ map1_val = map1_buf[c];
+ map2_val = map2_buf[c];
+ if (Rast_is_d_null_value(&map1_val) ||
+ Rast_is_d_null_value(&map2_val))
+ continue;
+
+ sumX += map1_val;
+ sumY += map2_val;
+ sumsqX += map1_val * map1_val;
+ sumsqY += map2_val * map2_val;
+ sumXY += map1_val * map2_val;
+ count++;
+ }
+ }
+
+ B = (sumXY - sumX * sumY / count) / (sumsqX - sumX * sumX / count);
+ R = (sumXY - sumX * sumY / count) /
+ sqrt((sumsqX - sumX * sumX / count) * (sumsqY - sumY * sumY / count));
+
+ meanX = sumX / count;
+ sumsqX = sumsqX / count;
+ varX = sumsqX - (meanX * meanX);
+ sdX = sqrt(varX);
+
+ meanY = sumY / count;
+ sumsqY = sumsqY / count;
+ varY = sumsqY - (meanY * meanY);
+ sdY = sqrt(varY);
+
+ A = meanY - B * meanX;
+ F = R * R / (1 - R * R / count - 2);
+
+ if (shell_style->answer) {
+ fprintf(stdout, "a=%f\n", A);
+ fprintf(stdout, "b=%f\n", B);
+ fprintf(stdout, "R=%f\n", R);
+ fprintf(stdout, "N=%ld\n", count);
+ fprintf(stdout, "F=%f\n", F);
+ fprintf(stdout, "meanX=%f\n", meanX);
+ fprintf(stdout, "sdX=%f\n", sdX);
+ fprintf(stdout, "meanY=%f\n", meanY);
+ fprintf(stdout, "sdY=%f\n", sdY);
+ }
+ else {
+ fprintf(stdout, "y = a + b*x\n");
+ fprintf(stdout, " a (Offset): %f\n", A);
+ fprintf(stdout, " b (Gain): %f\n", B);
+ fprintf(stdout, " R (sumXY - sumX*sumY/N): %f\n", R);
+ fprintf(stdout, " N (Number of elements): %ld\n", count);
+ fprintf(stdout, " F (F-test significance): %f\n", F);
+ fprintf(stdout, " meanX (Mean of map1): %f\n", meanX);
+ fprintf(stdout, " sdX (Standard deviation of map1): %f\n", sdX);
+ fprintf(stdout, " meanY (Mean of map2): %f\n", meanY);
+ fprintf(stdout, " sdY (Standard deviation of map1): %f", sdY);
+ }
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass/trunk/raster/r.regression.line/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/raster/r.regression.line/r.regression.line.html
===================================================================
--- grass/trunk/raster/r.regression.line/r.regression.line.html (rev 0)
+++ grass/trunk/raster/r.regression.line/r.regression.line.html 2010-09-24 15:24:55 UTC (rev 43677)
@@ -0,0 +1,49 @@
+<H2>DESCRIPTION</H2>
+
+<EM>r.regression.line</EM> Calculates linear regression from two raster maps,
+according to the formula y = a + b*x, where x and y represent raster maps.
+Optionally saves regression coefficients to an ASCII file.
+The result includes the following coefficients:
+offset/intercept (a) and gain/slope (b), residuals (R),
+number of elements (N), means (medX, medY), standard deviations
+(sdX, sdY), and the F test for testing the significance of the
+regression model as a whole (F).
+
+<H2>NOTES</H2>
+The flag <em>-s</em> select the slower method (applies to FP maps only)
+which writes out all pixel values individually to the temporary file.
+The results for offset/intercept (a) and gain/slope (b) are then
+identical to that obtained from R-stats's lm() function.
+
+<H2>EXAMPLE</H2>
+
+Comparison of the old and the new DEM in Spearfish:
+<div class="code"><pre>
+g.region rast=elevation.10m -p
+r.regression.line map1=elevation.dem map2=elevation.10m
+</pre></div>
+<p>
+
+Using the script style flag AND <em>eval</em> to make results
+available in the shell:
+<div class="code"><pre>
+g.region rast=elevation.10m -p
+eval `r.regression.line -g map1=elevation.dem map2=elevation.10m`
+echo $a
+479.615
+
+echo $b
+0.645631
+
+echo $R
+0.804441
+</pre></div>
+
+
+<H2>AUTHOR</H2>
+
+Dr. Agustin Lobo - alobo at ija.csic.es<BR>
+Updated to GRASS 5.7 Michael Barton, Arizona State University<BR>
+Script style output Markus Neteler
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass/trunk/raster/r.regression.line/r.regression.line.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list