1 /****************************************************************************
2 *
3 * MODULE: r.usler
4 * AUTHOR(S): Natalia Medvedeva - natmead@gmail.com
5 * Yann Chemin - yann.chemin@gmail.com
6 * PURPOSE: Calculates USLE R factor
7 * Rainfall Erosion index according to four methods
8 *
9 * COPYRIGHT: (C) 2002-2008, 2010 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * License (>=v2). Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 *****************************************************************************/
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <string.h>
20 #include <grass/gis.h>
21 #include <grass/raster.h>
22 #include <grass/glocale.h>
23
24 double elswaify_1985(double annual_pmm);
25 double morgan_1974(double annual_pmm);
26 double foster_1981(double annual_pmm);
27 double roose_1975(double annual_pmm);
28
main(int argc,char * argv[])29 int main(int argc, char *argv[])
30 {
31 int nrows, ncols;
32 int row, col;
33 char *nameflag; /*Switch for particular method */
34 struct GModule *module;
35 struct Option *input1, *input2, *output;
36 struct History history; /*metadata */
37 char *desc;
38
39 /************************************/
40 char *result; /*output raster name */
41 int infd_annual_pmm;
42 int outfd;
43 char *annual_pmm;
44
45 void *inrast_annual_pmm;
46 DCELL * outrast;
47
48 /************************************/
49 G_gisinit(argv[0]);
50
51 module = G_define_module();
52 G_add_keyword(_("raster"));
53 G_add_keyword(_("hydrology"));
54 G_add_keyword(_("rainfall"));
55 G_add_keyword(_("soil"));
56 G_add_keyword(_("erosion"));
57 module->description = _("Computes USLE R factor, Rainfall erosivity index.");
58
59 input2 = G_define_standard_option(G_OPT_R_INPUT);
60 input2->description = _("Name of annual precipitation raster map [mm/year]");
61
62 output = G_define_standard_option(G_OPT_R_OUTPUT);
63 output->description = _("Name for output USLE R raster map [MJ.mm/ha.hr.year]");
64
65 /* Define the different options */
66 input1 = G_define_option();
67 input1->key = "method";
68 input1->type = TYPE_STRING;
69 input1->required = YES;
70 input1->description = _("Name of USLE R equation");
71 input1->options = "roose, morgan, foster, elswaify";
72 desc = NULL;
73 G_asprintf(&desc,
74 "roose;%s;morgan;%s;foster;%s;elswaify;%s",
75 _("Roosle (1975)"),
76 _("Morgan (1974)"),
77 _("Foster (1981)"),
78 _("El-Swaify (1985)"));
79 input1->descriptions = desc;
80 input1->answer = "morgan";
81
82 /********************/
83 if (G_parser(argc, argv))
84 exit(EXIT_FAILURE);
85
86 nameflag = input1->answer;
87 annual_pmm = input2->answer;
88 result = output->answer;
89
90 /***************************************************/
91 infd_annual_pmm = Rast_open_old(annual_pmm, "");
92 inrast_annual_pmm = Rast_allocate_d_buf();
93
94 /***************************************************/
95 nrows = Rast_window_rows();
96 ncols = Rast_window_cols();
97 outrast = Rast_allocate_d_buf();
98
99 /* Create New raster files */
100 outfd = Rast_open_new(result, DCELL_TYPE);
101
102 /* Process pixels */
103 for (row = 0; row < nrows; row++)
104 {
105 DCELL d;
106 DCELL d_annual_pmm;
107 G_percent(row, nrows, 2);
108
109 /* read input map */
110 Rast_get_d_row(infd_annual_pmm, inrast_annual_pmm, row);
111
112 /*process the data */
113 for (col = 0; col < ncols; col++)
114 {
115 d_annual_pmm = ((DCELL *) inrast_annual_pmm)[col];
116 if (Rast_is_d_null_value(&d_annual_pmm))
117 Rast_set_d_null_value(&outrast[col], 1);
118 else
119 {
120 /*calculate morgan */
121 if (!strcmp(nameflag, "morgan"))
122 d = morgan_1974(d_annual_pmm);
123 /*calculate roose */
124 if (!strcmp(nameflag, "roose"))
125 d = roose_1975(d_annual_pmm);
126 /*calculate foster */
127 if (!strcmp(nameflag, "foster"))
128 d = foster_1981(d_annual_pmm);
129 /*calculate elswaify */
130 if (!strcmp(nameflag, "elswaify"))
131 d = elswaify_1985(d_annual_pmm);
132 outrast[col] = d ;
133 }
134 }
135 Rast_put_d_row(outfd, outrast);
136 }
137 G_free(inrast_annual_pmm);
138 Rast_close(infd_annual_pmm);
139 G_free(outrast);
140 Rast_close(outfd);
141
142 Rast_short_history(result, "raster", &history);
143 Rast_command_history(&history);
144 Rast_write_history(result, &history);
145
146 exit(EXIT_SUCCESS);
147 }
148