1 /****************************************************************************
2 *
3 * MODULE: r.uslek
4 * AUTHOR(S): Yann Chemin - yann.chemin@gmail.com
5 * PURPOSE: Transforms percentage of texture (sand/clay/silt)
6 * into USDA 1951 (p209) soil texture classes and then
7 * into USLE soil erodibility factor (K) as an output
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
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 #include <math.h>
22 #include <grass/gis.h>
23 #include <grass/raster.h>
24 #include <grass/glocale.h>
25
26 #define POLYGON_DIMENSION 20
27
28 double tex2usle_k(int texture, double om_in);
29 int prct2tex(double sand_input, double clay_input, double silt_input);
30 double point_in_triangle(double point_x, double point_y, double point_z,
31 double t1_x, double t1_y, double t1_z, double t2_x,
32 double t2_y, double t2_z, double t3_x, double t3_y,
33 double t3_z);
34
main(int argc,char * argv[])35 int main(int argc, char *argv[])
36 {
37 int nrows, ncols;
38 int row, col;
39 struct GModule *module;
40 struct Option *input1, *input2, *input3, *input4, *output1;
41 struct History history; /*metadata */
42
43 char *result; /*output raster name */
44 int infd_psand, infd_psilt, infd_pclay, infd_pomat;
45 int outfd;
46 char *psand, *psilt, *pclay, *pomat;
47 void *inrast_psand, *inrast_psilt, *inrast_pclay, *inrast_pomat;
48 DCELL *outrast;
49
50 /************************************/
51 G_gisinit(argv[0]);
52 module = G_define_module();
53 G_add_keyword(_("raster"));
54 G_add_keyword(_("hydrology"));
55 G_add_keyword(_("soil"));
56 G_add_keyword(_("erosion"));
57 module->description = _("Computes USLE Soil Erodibility Factor (K).");
58
59 /* Define the different options */
60 input1 = G_define_standard_option(G_OPT_R_INPUT);
61 input1->key = "psand";
62 input1->description = _("Name of soil sand fraction raster map [0.0-1.0]");
63
64 input2 = G_define_standard_option(G_OPT_R_INPUT);
65 input2->key = "pclay";
66 input2->description = _("Name of soil clay fraction raster map [0.0-1.0]");
67
68 input3 = G_define_standard_option(G_OPT_R_INPUT);
69 input3->key = "psilt";
70 input3->description = _("Name of soil silt fraction raster map [0.0-1.0]");
71
72 input4 = G_define_standard_option(G_OPT_R_INPUT);
73 input4->key = "pomat";
74 input4->description = _("Name of soil organic matter raster map [0.0-1.0]");
75
76 output1 = G_define_standard_option(G_OPT_R_OUTPUT);
77 output1->description = _("Name for output USLE K factor raster map [t.ha.hr/ha.MJ.mm]");
78
79 /********************/
80 if (G_parser(argc, argv))
81 exit(EXIT_FAILURE);
82
83 psand = input1->answer;
84 pclay = input2->answer;
85 psilt = input3->answer;
86 pomat = input4->answer;
87 result = output1->answer;
88
89 /***************************************************/
90 infd_psand = Rast_open_old(psand, "");
91 inrast_psand = Rast_allocate_d_buf();
92
93 infd_psilt = Rast_open_old(psilt, "");
94 inrast_psilt = Rast_allocate_d_buf();
95
96 infd_pclay = Rast_open_old(pclay, "");
97 inrast_pclay = Rast_allocate_d_buf();
98
99 infd_pomat = Rast_open_old(pomat, "");
100 inrast_pomat = Rast_allocate_d_buf();
101 /***************************************************/
102 nrows = Rast_window_rows();
103 ncols = Rast_window_cols();
104 outrast = Rast_allocate_d_buf();
105
106 /* Create New raster files */
107 outfd = Rast_open_new(result, DCELL_TYPE);
108
109 /* Process pixels */
110 for (row = 0; row < nrows; row++)
111 {
112 DCELL d;
113 DCELL d_sand;
114 DCELL d_clay;
115 DCELL d_silt;
116 DCELL d_om;
117 G_percent(row, nrows, 2);
118
119 /* read soil input maps */
120 Rast_get_d_row(infd_psand, inrast_psand, row);
121 Rast_get_d_row(infd_psilt, inrast_psilt, row);
122 Rast_get_d_row(infd_pclay, inrast_pclay, row);
123 Rast_get_d_row(infd_pomat, inrast_pomat, row);
124
125 /*process the data */
126 for (col = 0; col < ncols; col++)
127 {
128 d_sand = ((DCELL *) inrast_psand)[col];
129 d_silt = ((DCELL *) inrast_psilt)[col];
130 d_clay = ((DCELL *) inrast_pclay)[col];
131 d_om = ((DCELL *) inrast_pomat)[col];
132 if (Rast_is_d_null_value(&d_sand) ||
133 Rast_is_d_null_value(&d_clay) ||
134 Rast_is_d_null_value(&d_silt))
135 Rast_set_d_null_value(&outrast[col], 1);
136 else {
137 /***************************************/
138 /* In case some map input not standard */
139 if (fabs(d_sand + d_clay + d_silt - 1.0) > GRASS_EPSILON )
140 Rast_set_d_null_value(&outrast[col], 1);
141 /* if OM==NULL then make it 0.0 */
142 else {
143 if (Rast_is_d_null_value(&d_om))
144 d_om = 0.0;
145 /************************************/
146 /* convert to usle_k */
147 d = (double)prct2tex(d_sand, d_clay, d_silt);
148 if(d > 11.0) /* 11 is highest class */
149 Rast_set_d_null_value(&outrast[col], 1);
150 else
151 outrast[col] = tex2usle_k((int)d, d_om);
152 }
153 }
154 }
155 Rast_put_d_row(outfd, outrast);
156 }
157 G_free(inrast_psand);
158 G_free(inrast_psilt);
159 G_free(inrast_pclay);
160 G_free(inrast_pomat);
161 Rast_close(infd_psand);
162 Rast_close(infd_psilt);
163 Rast_close(infd_pclay);
164 Rast_close(infd_pomat);
165 G_free(outrast);
166 Rast_close(outfd);
167
168 Rast_short_history(result, "raster", &history);
169 Rast_command_history(&history);
170 Rast_write_history(result, &history);
171
172 exit(EXIT_SUCCESS);
173 }
174