← Back to branch summary

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

« back to all changes in this revision

Viewing changes to imagery/i.topo.corr/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.topo.corr/main.c
19
19
#include <stdlib.h>
20
20
#include <string.h>
21
21
#include <grass/gis.h>
 
22
#include <grass/raster.h>
22
23
#include <grass/glocale.h>
23
24
 
24
25
#include "local_proto.h"
25
26
 
26
 
int full_open_old(Gfile * gf, char *fname)
27
 
{
28
 
    gf->fd = -1;
29
 
 
30
 
    snprintf(gf->name, 127, "%s", fname);
31
 
 
32
 
    gf->mapset = G_find_cell2(gf->name, "");
33
 
    if (gf->mapset == NULL)
34
 
        G_fatal_error(_("Raster map <%s> not found"), gf->name);
35
 
 
36
 
    gf->fd = G_open_cell_old(gf->name, gf->mapset);
37
 
    if (gf->fd < 0)
38
 
        G_fatal_error(_("Unable to open raster map <%s@%s>"), gf->name,
39
 
                      gf->mapset);
40
 
 
41
 
    gf->type = G_raster_map_type(gf->name, gf->mapset);
42
 
 
43
 
    return gf->fd;
44
 
}
45
 
 
46
 
int full_open_new(Gfile * gf, char *fname, RASTER_MAP_TYPE ftype)
47
 
{
48
 
    gf->fd = -1;
49
 
 
50
 
    snprintf(gf->name, 127, "%s", fname);
51
 
    if (G_legal_filename(gf->name) < 0)
52
 
        G_fatal_error(_("<%s> is an illegal name"), gf->name);
53
 
 
54
 
    gf->type = ftype;
55
 
    gf->mapset = G_mapset();
56
 
    gf->fd = G_open_raster_new(gf->name, gf->type);
57
 
    if (gf->fd < 0)
58
 
        G_fatal_error(_("Unable to create raster map <%s>"), gf->name);
59
 
    
60
 
    return gf->fd;
61
 
}
62
 
 
63
27
int main(int argc, char *argv[])
64
28
{
65
29
    struct History history;
66
30
    struct GModule *module;
67
31
    struct Cell_head hd_band, hd_dem, window;
68
32
 
69
 
    char bufname[128];          /* TODO: use GNAME_MAX? */
70
 
 
71
33
    int i;
72
34
    struct Option *base, *output, *input, *zeni, *azim, *metho;
73
 
    struct Flag *ilum;
 
35
    struct Flag *ilum, *scl;
74
36
 
75
37
    Gfile dem, out, band;
76
38
    double zenith, azimuth;
77
39
    int method = COSINE;
 
40
    int do_scale;
78
41
 
79
42
    /* initialize GIS environment */
80
43
    G_gisinit(argv[0]);
81
44
 
82
45
    /* initialize module */
83
46
    module = G_define_module();
 
47
    G_add_keyword(_("imagery"));
 
48
    G_add_keyword(_("terrain"));
 
49
    G_add_keyword(_("topographic correction"));
84
50
    module->description = _("Computes topographic correction of reflectance.");
85
 
    module->keywords =
86
 
        _("imagery, terrain, topographic correction");
87
 
    
 
51
 
88
52
    /* It defines the different parameters */
89
53
 
90
54
    input = G_define_standard_option(G_OPT_R_INPUTS);
125
89
    ilum->key = 'i';
126
90
    ilum->description = _("Output sun illumination terrain model");
127
91
    
 
92
    scl = G_define_flag();
 
93
    scl->key = 's';
 
94
    scl->description = _("Scale output to input and copy color rules");
 
95
    
128
96
    if (G_parser(argc, argv))
129
97
        exit(EXIT_FAILURE);
130
98
 
132
100
        G_fatal_error(_("Solar azimuth is necessary to calculate illumination terrain model"));
133
101
 
134
102
    if (!ilum->answer && input->answer == NULL)
135
 
        G_fatal_error(_("Reflectance maps are necessary to make topographic correction"));
 
103
        G_fatal_error
 
104
            (_("Reflectance maps are necessary to make topographic correction"));
136
105
 
137
106
    zenith = atof(zeni->answer);
 
107
    out.type = DCELL_TYPE;
 
108
    do_scale = scl->answer;
138
109
 
139
110
    /* Evaluate only cos_i raster file */
140
111
    /* i.topo.corr -i out=cosi.on07 base=SRTM_v2 zenith=33.3631 azimuth=59.8897 */
141
112
    if (ilum->answer) {
 
113
        Rast_get_window(&window);
142
114
        azimuth = atof(azim->answer);
143
115
        /* Warning: make buffers and output after set window */
144
 
        full_open_old(&dem, base->answer);
 
116
        strcpy(dem.name, base->answer);
145
117
        /* Set window to DEM file */
146
 
        G_get_window(&window);
147
 
        G_get_cellhd(dem.name, dem.mapset, &hd_dem);
148
 
        G_align_window(&window, &hd_dem);
149
 
        G_set_window(&window);
 
118
        Rast_get_window(&window);
 
119
        Rast_get_cellhd(dem.name, "", &hd_dem);
 
120
        Rast_align_window(&window, &hd_dem);
 
121
        dem.fd = Rast_open_old(dem.name, "");
 
122
        dem.type = Rast_get_map_type(dem.fd);
150
123
        /* Open and buffer of the output file */
151
 
        full_open_new(&out, output->answer, DCELL_TYPE);
152
 
        out.rast = G_allocate_raster_buf(out.type);
 
124
        strcpy(out.name, output->answer);
 
125
        out.fd = Rast_open_new(output->answer, DCELL_TYPE);
 
126
        out.rast = Rast_allocate_buf(out.type);
153
127
        /* Open and buffer of the elevation file */
154
 
        if (dem.type == CELL_TYPE) {
155
 
            dem.rast = G_allocate_raster_buf(CELL_TYPE);
156
 
            eval_c_cosi(&out, &dem, zenith, azimuth);
157
 
        }
158
 
        else if (dem.type == FCELL_TYPE) {
159
 
            dem.rast = G_allocate_raster_buf(FCELL_TYPE);
160
 
            eval_f_cosi(&out, &dem, zenith, azimuth);
161
 
        }
162
 
        else if (dem.type == DCELL_TYPE) {
163
 
            dem.rast = G_allocate_raster_buf(DCELL_TYPE);
164
 
            eval_d_cosi(&out, &dem, zenith, azimuth);
165
 
        }
166
 
        else {
167
 
            G_fatal_error(_("Elevation raster map of unknown type"));
168
 
        }
 
128
        dem.rast = Rast_allocate_buf(dem.type);
 
129
        eval_cosi(&out, &dem, zenith, azimuth);
169
130
        /* Close files, buffers, and write history */
170
131
        G_free(dem.rast);
171
 
        G_close_cell(dem.fd);
 
132
        Rast_close(dem.fd);
172
133
        G_free(out.rast);
173
 
        G_close_cell(out.fd);
174
 
        G_short_history(out.name, "raster", &history);
175
 
        G_command_history(&history);
176
 
        G_write_history(out.name, &history);
 
134
        Rast_close(out.fd);
 
135
        Rast_short_history(out.name, "raster", &history);
 
136
        Rast_command_history(&history);
 
137
        Rast_write_history(out.name, &history);
177
138
    }
178
139
    /* Evaluate topographic correction for all bands */
179
140
    /* i.topo.corr input=on07.toar.1 out=tcor base=cosi.on07 zenith=33.3631 method=c-factor */
196
157
        else
197
158
            G_fatal_error(_("Invalid method: %s"), metho->answer);
198
159
 
199
 
        full_open_old(&dem, base->answer);
 
160
        dem.fd = Rast_open_old(base->answer, "");
 
161
        dem.type = Rast_get_map_type(dem.fd);
 
162
        Rast_close(dem.fd);
200
163
        if (dem.type == CELL_TYPE)
201
164
            G_fatal_error(_("Illumination model is of CELL type"));
202
165
 
203
166
        for (i = 0; input->answers[i] != NULL; i++) {
204
 
            G_message("Band %s: ", input->answers[i]);
 
167
            G_message(_("Band %s: "), input->answers[i]);
205
168
            /* Abre fichero de bandas y el de salida */
206
 
            full_open_old(&band, input->answers[i]);
 
169
            strcpy(band.name, input->answers[i]);
 
170
            Rast_get_cellhd(band.name, "", &hd_band);
 
171
            Rast_set_window(&hd_band);  /* Antes de out_open y allocate para mismo size */
 
172
            band.fd = Rast_open_old(band.name, "");
 
173
            band.type = Rast_get_map_type(band.fd);
207
174
            if (band.type != DCELL_TYPE) {
208
 
                G_warning(_("Reflectance of raster map <%s> is not of DCELL type - ignored"),
 
175
                G_warning(_("Reflectance of <%s> is not of DCELL type - ignored."),
209
176
                          input->answers[i]);
210
 
                G_close_cell(band.fd);
 
177
                Rast_close(band.fd);
211
178
                continue;
212
179
            }
213
 
            G_get_cellhd(band.name, band.mapset, &hd_band);
214
 
            G_set_window(&hd_band);     /* Antes de out_open y allocate para mismo tama�o */
215
180
            /* ----- */
216
 
            snprintf(bufname, 127, "%s.%s", output->answer,
 
181
            dem.fd = Rast_open_old(base->answer, "");
 
182
            snprintf(out.name, GNAME_MAX - 1, "%s.%s", output->answer,
217
183
                     input->answers[i]);
218
 
            full_open_new(&out, bufname, DCELL_TYPE);
219
 
            out.rast = G_allocate_raster_buf(out.type);
220
 
            band.rast = G_allocate_raster_buf(band.type);
221
 
            dem.rast = G_allocate_raster_buf(dem.type);
 
184
            out.fd = Rast_open_new(out.name, DCELL_TYPE);
 
185
            out.rast = Rast_allocate_buf(out.type);
 
186
            band.rast = Rast_allocate_buf(band.type);
 
187
            dem.rast = Rast_allocate_buf(dem.type);
222
188
            /* ----- */
223
 
            eval_tcor(method, &out, &dem, &band, zenith);
 
189
            eval_tcor(method, &out, &dem, &band, zenith, do_scale);
224
190
            /* ----- */
225
191
            G_free(dem.rast);
 
192
            Rast_close(dem.fd);
226
193
            G_free(band.rast);
227
 
            G_close_cell(band.fd);
 
194
            Rast_close(band.fd);
228
195
            G_free(out.rast);
229
 
            G_close_cell(out.fd);
230
 
            G_short_history(out.name, "raster", &history);
231
 
            G_command_history(&history);
232
 
            G_write_history(out.name, &history);
233
 
 
234
 
            char command[300];
235
 
 
236
 
            /* TODO: better avoid system() */
237
 
            sprintf(command, "r.colors map=%s color=grey", out.name);
238
 
            system(command);
239
 
 
240
 
/* new but not functional:
 
196
            Rast_close(out.fd);
 
197
            Rast_short_history(out.name, "raster", &history);
 
198
            Rast_command_history(&history);
 
199
            Rast_write_history(out.name, &history);
 
200
 
241
201
            {
242
202
                struct FPRange range;
243
203
                DCELL min, max;
244
204
                struct Colors grey;
245
 
                G_read_fp_range(out.name, G_mapset(), &range);
246
 
                G_get_fp_range_min_max(&range, &min, &max);
247
 
                G_make_grey_scale_colors(&grey, min, max);
248
 
                G_write_colors(out.name, G_mapset(), &grey);
 
205
                int make_colors = 1;
 
206
 
 
207
                if (do_scale) {
 
208
                    if (Rast_read_colors(band.name, "", &grey) >= 0)
 
209
                        make_colors = 0;
 
210
                }
 
211
                
 
212
                if (make_colors) {
 
213
                    Rast_read_fp_range(out.name, G_mapset(), &range);
 
214
                    Rast_get_fp_range_min_max(&range, &min, &max);
 
215
                    Rast_make_grey_scale_colors(&grey, min, max);
 
216
                }
 
217
                Rast_write_colors(out.name, G_mapset(), &grey);
249
218
            }
250
 
*/
251
 
 
252
219
        }
253
 
        G_close_cell(dem.fd);
254
220
    }
255
221
 
256
222
    exit(EXIT_SUCCESS);