← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to imagery/i.aster.toar/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
imagery/i.aster.toar/main.c
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       i.aster.toar
 
5
 * AUTHOR(S):    Yann Chemin - yann.chemin@gmail.com
 
6
 * PURPOSE:      Calculate TOA Reflectance for Aster from DN.
 
7
 *               Input 9 bands (VNIR and SWIR).
 
8
 *
 
9
 * COPYRIGHT:    (C) 2002-2010 by the GRASS Development Team
 
10
 *
 
11
 *               This program is free software under the GNU Lesser General Public
 
12
 *               License. Read the file COPYING that comes with GRASS for details.
 
13
 *
 
14
 *****************************************************************************/
 
15
 
 
16
#include <stdio.h>
 
17
#include <stdlib.h>
 
18
#include <string.h>
 
19
#include <grass/gis.h>
 
20
#include <grass/raster.h>
 
21
#include <grass/glocale.h>
 
22
 
 
23
#define MAXFILES 15
 
24
 
 
25
/* DN to radiance conversion factors */
 
26
double gain_aster(int band_number, int gain_code);
 
27
 
 
28
    /*Gain Code */
 
29
    /*0 - High (Not Applicable for band 10-14: TIR) */
 
30
    /*1 - Normal */
 
31
    /*2 - Low 1(Not Applicable for band 10-14: TIR) */
 
32
    /*3 - Low 2(Not Applicable for Band 1-3N/B & 10-14) */
 
33
 
 
34
/*sun exo-atmospheric irradiance */
 
35
#define KEXO1 1828.0
 
36
#define KEXO2 1559.0
 
37
#define KEXO3 1045.0
 
38
#define KEXO4 226.73
 
39
#define KEXO5 86.50
 
40
#define KEXO6 81.99
 
41
#define KEXO7 74.72
 
42
#define KEXO8 66.41
 
43
#define KEXO9 59.83
 
44
 
 
45
#define PI M_PI
 
46
 
 
47
double rad2ref_aster(double radiance, double doy, double sun_elevation,
 
48
                     double k_exo);
 
49
 
 
50
int main(int argc, char *argv[])
 
51
{
 
52
    struct Cell_head cellhd;    /*region+header info */
 
53
    char *mapset;               /*mapset name */
 
54
    int nrows, ncols;
 
55
    int row, col;
 
56
    struct GModule *module;
 
57
    struct Option *input, *output;
 
58
    struct Option *input1, *input2;
 
59
    struct Flag *flag0, *flag1, *flag2;
 
60
    struct Flag *flag3, *flag4, *flag5;
 
61
    struct History history;     /*metadata */
 
62
 
 
63
    /************************************/
 
64
    char *name;                 /*input raster name */
 
65
    char *result;               /*output raster name */
 
66
 
 
67
    /*Prepare new names for output files */
 
68
    char result0[GNAME_MAX], result1[GNAME_MAX];
 
69
    char result2[GNAME_MAX], result3[GNAME_MAX];
 
70
    char result4[GNAME_MAX], result5[GNAME_MAX];
 
71
    char result6[GNAME_MAX], result7[GNAME_MAX];
 
72
    char result8[GNAME_MAX], result9[GNAME_MAX];
 
73
    char result10[GNAME_MAX], result11[GNAME_MAX];
 
74
    char result12[GNAME_MAX], result13[GNAME_MAX];
 
75
    char result14[GNAME_MAX];
 
76
 
 
77
    /*File Descriptors */
 
78
    int infd[MAXFILES];
 
79
    int outfd[MAXFILES];
 
80
    char **names, **ptr;
 
81
    /* For some strange reason infd[0] cannot be used later */
 
82
    /* So nfiles is initialized with nfiles = 1 */
 
83
    int nfiles = 1;
 
84
    int i = 0, j = 0;
 
85
    int radiance = 0;
 
86
    void *inrast[MAXFILES];
 
87
    DCELL *outrast[MAXFILES];
 
88
    RASTER_MAP_TYPE in_data_type[MAXFILES];
 
89
    RASTER_MAP_TYPE out_data_type = DCELL_TYPE; /* 0=numbers  1=text */
 
90
    double gain[MAXFILES], offset[MAXFILES];
 
91
    double kexo[MAXFILES];
 
92
    double doy, sun_elevation;
 
93
 
 
94
    /************************************/
 
95
    G_gisinit(argv[0]);
 
96
    module = G_define_module();
 
97
    G_add_keyword(_("imagery"));
 
98
    G_add_keyword(_("radiometric conversion"));
 
99
    G_add_keyword(_("Terra-ASTER"));
 
100
    G_add_keyword(_("radiance"));
 
101
    G_add_keyword(_("reflectance"));
 
102
    G_add_keyword(_("brightness temperature"));
 
103
    module->description =
 
104
        _("Calculates Top of Atmosphere Radiance/Reflectance/Brightness Temperature from ASTER DN.\n");
 
105
 
 
106
    /* Define the different options */
 
107
    input = G_define_standard_option(G_OPT_R_INPUTS);
 
108
    input->description = _("Names of ASTER DN layers (15 layers)");
 
109
 
 
110
    input1 = G_define_option();
 
111
    input1->key = "dayofyear";
 
112
    input1->type = TYPE_DOUBLE;
 
113
    input1->required = YES;
 
114
    input1->gisprompt = "value";
 
115
    input1->description = _("Day of Year of satellite overpass [0-366]");
 
116
 
 
117
    input2 = G_define_option();
 
118
    input2->key = "sun_elevation";
 
119
    input2->type = TYPE_DOUBLE;
 
120
    input2->required = YES;
 
121
    input2->gisprompt = "value";
 
122
    input2->description = _("Sun elevation angle (degrees, < 90.0)");
 
123
 
 
124
    output = G_define_standard_option(G_OPT_R_OUTPUT);
 
125
    output->description = _("Base name of the output layers (will add .x)");
 
126
 
 
127
    /* Define the different flags */
 
128
    flag0 = G_define_flag();
 
129
    flag0->key = 'r';
 
130
    flag0->description = _("Output is radiance (W/m2)");
 
131
 
 
132
    flag1 = G_define_flag();
 
133
    flag1->key = 'a';
 
134
    flag1->description = _("VNIR is High Gain");
 
135
 
 
136
    flag2 = G_define_flag();
 
137
    flag2->key = 'b';
 
138
    flag2->description = _("SWIR is High Gain");
 
139
 
 
140
    flag3 = G_define_flag();
 
141
    flag3->key = 'c';
 
142
    flag3->description = _("VNIR is Low Gain 1");
 
143
 
 
144
    flag4 = G_define_flag();
 
145
    flag4->key = 'd';
 
146
    flag4->description = _("SWIR is Low Gain 1");
 
147
 
 
148
    flag5 = G_define_flag();
 
149
    flag5->key = 'e';
 
150
    flag5->description = _("SWIR is Low Gain 2");
 
151
 
 
152
    /********************/
 
153
    if (G_parser(argc, argv))
 
154
        exit(EXIT_FAILURE);
 
155
 
 
156
    names = input->answers;
 
157
    ptr = input->answers;
 
158
    doy = atof(input1->answer);
 
159
    sun_elevation = atof(input2->answer);
 
160
    result = output->answer;
 
161
 
 
162
    radiance = (flag0->answer);
 
163
 
 
164
    /********************/
 
165
    /*Prepare the output file names */
 
166
 
 
167
    /********************/
 
168
    sprintf(result0,"%s%s", result, ".1");
 
169
    sprintf(result1,"%s%s", result, ".2");
 
170
    sprintf(result2,"%s%s", result, ".3N");
 
171
    sprintf(result3,"%s%s", result, ".3B");
 
172
    sprintf(result4,"%s%s", result, ".4");
 
173
    sprintf(result5,"%s%s", result, ".5");
 
174
    sprintf(result6,"%s%s", result, ".6");
 
175
    sprintf(result7,"%s%s", result, ".7");
 
176
    sprintf(result8,"%s%s", result, ".8");
 
177
    sprintf(result9,"%s%s", result, ".9");
 
178
    sprintf(result10,"%s%s", result, ".10");
 
179
    sprintf(result11,"%s%s", result, ".11");
 
180
    sprintf(result12,"%s%s", result, ".12");
 
181
    sprintf(result13,"%s%s", result, ".13");
 
182
    sprintf(result14,"%s%s", result, ".14");
 
183
    /********************/
 
184
    /*Prepare radiance boundaries */
 
185
 
 
186
    /********************/
 
187
    int gain_code = 1;
 
188
 
 
189
    for (i = 0; i < MAXFILES; i++) {
 
190
        /*0 - High (Not Applicable for band 10-14: TIR) */
 
191
        /*1 - Normal */
 
192
        /*2 - Low 1(Not Applicable for band 10-14: TIR) */
 
193
        /*3 - Low 2(Not Applicable for Band 1-3N/B & 10-14) */
 
194
        if (flag1->answer && i <= 3)
 
195
            gain_code = 0;
 
196
        if (flag2->answer && i >= 4 && i <= 9)
 
197
            gain_code = 0;
 
198
        if (flag3->answer && i <= 3)
 
199
            gain_code = 2;
 
200
        if (flag4->answer && i >= 4 && i <= 9)
 
201
            gain_code = 2;
 
202
        if (flag5->answer && i >= 4 && i <= 9)
 
203
            gain_code = 3;
 
204
        gain[i] = gain_aster(i, gain_code);
 
205
        /* Reset to NORMAL GAIN */
 
206
        gain_code = 1;
 
207
    }
 
208
 
 
209
    /********************/
 
210
    /*Prepare sun exo-atm irradiance */
 
211
 
 
212
    /********************/
 
213
    kexo[0] = KEXO1;
 
214
    kexo[1] = KEXO2;
 
215
    kexo[2] = KEXO3;
 
216
    kexo[3] = KEXO3;
 
217
    kexo[4] = KEXO4;
 
218
    kexo[5] = KEXO5;
 
219
    kexo[6] = KEXO6;
 
220
    kexo[7] = KEXO7;
 
221
    kexo[8] = KEXO8;
 
222
    kexo[9] = KEXO9;
 
223
 
 
224
    /********************/
 
225
 
 
226
    /********************/
 
227
    for (; *ptr != NULL; ptr++) {
 
228
        if (nfiles > MAXFILES)
 
229
            G_fatal_error(_("Too many input maps. Only %d allowed."),
 
230
                          MAXFILES);
 
231
        name = *ptr;
 
232
        /* Allocate input buffer */
 
233
        in_data_type[nfiles-1] = Rast_map_type(name, "");
 
234
        /* For some strange reason infd[0] cannot be used later */
 
235
        /* So nfiles is initialized with nfiles = 1 */
 
236
        infd[nfiles] = Rast_open_old(name, "");
 
237
 
 
238
        Rast_get_cellhd(name, "", &cellhd);
 
239
        inrast[nfiles-1] = Rast_allocate_buf(in_data_type[nfiles-1]);
 
240
        nfiles++;
 
241
    }
 
242
    nfiles--;
 
243
    if (nfiles < MAXFILES) 
 
244
        G_fatal_error(_("The input band number should be 15"));
 
245
 
 
246
    /***************************************************/
 
247
    /* Allocate output buffer, use input map data_type */
 
248
    nrows = Rast_window_rows();
 
249
    ncols = Rast_window_cols();
 
250
    out_data_type = DCELL_TYPE;
 
251
    for (i = 0; i < MAXFILES; i++)
 
252
        outrast[i] = Rast_allocate_buf(out_data_type);
 
253
 
 
254
    outfd[1] = Rast_open_new(result0, 1);
 
255
    outfd[2] = Rast_open_new(result1, 1);
 
256
    outfd[3] = Rast_open_new(result2, 1);
 
257
    outfd[4] = Rast_open_new(result3, 1);
 
258
    outfd[5] = Rast_open_new(result4, 1);
 
259
    outfd[6] = Rast_open_new(result5, 1);
 
260
    outfd[7] = Rast_open_new(result6, 1);
 
261
    outfd[8] = Rast_open_new(result7, 1);
 
262
    outfd[9] = Rast_open_new(result8, 1);
 
263
    outfd[10] = Rast_open_new(result9, 1);
 
264
    outfd[11] = Rast_open_new(result10, 1);
 
265
    outfd[12] = Rast_open_new(result11, 1);
 
266
    outfd[13] = Rast_open_new(result12, 1);
 
267
    outfd[14] = Rast_open_new(result13, 1);
 
268
    outfd[15] = Rast_open_new(result14, 1);
 
269
    /* Process pixels */
 
270
 
 
271
    DCELL dout[MAXFILES];
 
272
    DCELL d[MAXFILES];
 
273
 
 
274
    for (row = 0; row < nrows; row++) {
 
275
        G_percent(row, nrows, 2);
 
276
        /* read input map */
 
277
        for (i = 1; i <= MAXFILES; i++)
 
278
            Rast_get_row(infd[i], inrast[i-1], row, in_data_type[i-1]);
 
279
 
 
280
        /*process the data */
 
281
        for (col = 0; col < ncols; col++) {
 
282
            for (i = 0; i < MAXFILES; i++) {
 
283
                switch (in_data_type[i]) {
 
284
                case CELL_TYPE:
 
285
                    d[i] = (double)((CELL *) inrast[i])[col];
 
286
                    break;
 
287
                case FCELL_TYPE:
 
288
                    d[i] = (double)((FCELL *) inrast[i])[col];
 
289
                    break;
 
290
                case DCELL_TYPE:
 
291
                    d[i] = (double)((DCELL *) inrast[i])[col];
 
292
                    break;
 
293
                }
 
294
                /* if radiance mode or Thermal band */
 
295
                if (radiance || i >= 10) {
 
296
                    dout[i] = gain[i] * (d[i] - 1.0);
 
297
                }
 
298
                /* if reflectance default mode and Not Thermal Band */
 
299
                else {
 
300
                    dout[i] = gain[i] * (d[i] - 1.0);
 
301
                    dout[i] =
 
302
                        rad2ref_aster(dout[i], doy, sun_elevation, kexo[i]);
 
303
                }
 
304
                outrast[i][col] = dout[i];
 
305
            }
 
306
        }
 
307
        for (i = 1; i <= MAXFILES; i++)
 
308
            Rast_put_row(outfd[i], outrast[i-1], out_data_type);
 
309
    }
 
310
    for (i = 1; i <= MAXFILES; i++) {
 
311
        G_free(inrast[i-1]);
 
312
        Rast_close(infd[i]);
 
313
        G_free(outrast[i-1]);
 
314
        Rast_close(outfd[i]);
 
315
    }
 
316
    exit(EXIT_SUCCESS);
 
317
}