← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to raster/r.sunhours/main.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.sunhours/main.c
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       r.sunhours
 
5
 * AUTHOR(S):    Markus Metz
 
6
 * PURPOSE:      Calculates solar azimuth and angle, and 
 
7
 *               sunshine hours (also called daytime period)
 
8
 *               Uses NREL SOLPOS
 
9
 * COPYRIGHT:    (C) 2010-2013 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
/* TODO: always use solpos if time is Greenwich standard time */
 
18
 
 
19
#define MAIN
 
20
 
 
21
#include <stdio.h>
 
22
#include <stdlib.h>
 
23
#include <string.h>
 
24
#include <math.h>
 
25
#include <grass/gis.h>
 
26
#include <grass/raster.h>
 
27
#include <grass/gprojects.h>
 
28
#include <grass/glocale.h>
 
29
#include "solpos00.h"
 
30
 
 
31
void set_solpos_time(struct posdata *, int, int, int, int, int, int);
 
32
void set_solpos_longitude(struct posdata *, double );
 
33
int roundoff(double *);
 
34
 
 
35
int main(int argc, char *argv[])
 
36
{
 
37
    struct GModule *module;
 
38
    struct
 
39
    {
 
40
        struct Option *elev, *azimuth, *sunhours, *year,
 
41
            *month, *day, *hour, *minutes, *seconds;
 
42
        struct Flag *lst_time, *no_solpos;
 
43
    } parm;
 
44
    struct Cell_head window;
 
45
    FCELL *elevbuf, *azimuthbuf, *sunhourbuf;
 
46
    struct History hist;
 
47
 
 
48
    /* projection information of input map */
 
49
    struct Key_Value *in_proj_info, *in_unit_info;
 
50
    struct pj_info iproj;       /* input map proj parameters  */
 
51
    struct pj_info oproj;       /* output map proj parameters  */
 
52
 
 
53
    char *elev_name, *azimuth_name, *sunhour_name;
 
54
    int elev_fd, azimuth_fd, sunhour_fd;
 
55
    double ha, ha_cos, s_gamma, s_elevation, s_azimuth;
 
56
    double s_declination, sd_sin, sd_cos;
 
57
    double se_sin, sa_cos;
 
58
    double east, east_ll, north, north_ll;
 
59
    double north_gc, north_gc_sin, north_gc_cos;  /* geocentric latitude */
 
60
    double ba2;
 
61
    int year, month, day, hour, minutes, seconds;
 
62
    int doy;                                    /* day of year */
 
63
    int row, col, nrows, ncols;
 
64
    int do_reproj = 0;
 
65
    int lst_time = 1;
 
66
    int use_solpos = 0;
 
67
    struct posdata pd;
 
68
 
 
69
    G_gisinit(argv[0]);
 
70
 
 
71
    module = G_define_module();
 
72
    G_add_keyword(_("raster"));
 
73
    G_add_keyword(_("solar"));
 
74
    module->label = _("Calculates solar elevation, solar azimuth, and sun hours.");
 
75
    module->description = _("Solar elevation: the angle between the direction of the geometric center "
 
76
                            "of the sun's apparent disk and the (idealized) horizon. "
 
77
                            "Solar azimuth: the angle from due north in clockwise direction.");
 
78
    
 
79
    parm.elev = G_define_standard_option(G_OPT_R_OUTPUT);
 
80
    parm.elev->key = "elevation";
 
81
    parm.elev->label = _("Output raster map with solar elevation angle");
 
82
    parm.elev->required = NO;
 
83
 
 
84
    parm.azimuth = G_define_standard_option(G_OPT_R_OUTPUT);
 
85
    parm.azimuth->key = "azimuth";
 
86
    parm.azimuth->label = _("Output raster map with solar azimuth angle");
 
87
    parm.azimuth->required = NO;
 
88
    
 
89
    parm.sunhours = G_define_standard_option(G_OPT_R_OUTPUT);
 
90
    parm.sunhours->key = "sunhour";
 
91
    parm.sunhours->label = _("Output raster map with sunshine hours");
 
92
    parm.sunhours->description = _("Sunshine hours require SOLPOS use and Greenwich standard time");
 
93
    parm.sunhours->required = NO;
 
94
 
 
95
    parm.year = G_define_option();
 
96
    parm.year->key = "year";
 
97
    parm.year->type = TYPE_INTEGER;
 
98
    parm.year->required = YES;
 
99
    parm.year->description = _("Year");
 
100
    parm.year->options = "1950-2050";
 
101
    parm.year->guisection = _("Time");
 
102
 
 
103
    parm.month = G_define_option();
 
104
    parm.month->key = "month";
 
105
    parm.month->type = TYPE_INTEGER;
 
106
    parm.month->required = NO;
 
107
    parm.month->label = _("Month");
 
108
    parm.month->description = _("If not given, day is interpreted as day of the year");
 
109
    parm.month->options = "1-12";
 
110
    parm.month->guisection = _("Time");
 
111
 
 
112
    parm.day = G_define_option();
 
113
    parm.day->key = "day";
 
114
    parm.day->type = TYPE_INTEGER;
 
115
    parm.day->required = YES;
 
116
    parm.day->description = _("Day");
 
117
    parm.day->options = "1-366";
 
118
    parm.day->guisection = _("Time");
 
119
 
 
120
    parm.hour = G_define_option();
 
121
    parm.hour->key = "hour";
 
122
    parm.hour->type = TYPE_INTEGER;
 
123
    parm.hour->required = NO;
 
124
    parm.hour->description = _("Hour");
 
125
    parm.hour->options = "0-24";
 
126
    parm.hour->answer = "12";
 
127
    parm.hour->guisection = _("Time");
 
128
 
 
129
    parm.minutes = G_define_option();
 
130
    parm.minutes->key = "minute";
 
131
    parm.minutes->type = TYPE_INTEGER;
 
132
    parm.minutes->required = NO;
 
133
    parm.minutes->description = _("Minutes");
 
134
    parm.minutes->options = "0-60";
 
135
    parm.minutes->answer = "0";
 
136
    parm.minutes->guisection = _("Time");
 
137
 
 
138
    parm.seconds = G_define_option();
 
139
    parm.seconds->key = "second";
 
140
    parm.seconds->type = TYPE_INTEGER;
 
141
    parm.seconds->required = NO;
 
142
    parm.seconds->description = _("Seconds");
 
143
    parm.seconds->options = "0-60";
 
144
    parm.seconds->answer = "0";
 
145
    parm.seconds->guisection = _("Time");
 
146
 
 
147
    parm.lst_time = G_define_flag();
 
148
    parm.lst_time->key = 't';
 
149
    parm.lst_time->description = _("Time is local sidereal time, not Greenwich standard time");
 
150
 
 
151
    parm.no_solpos = G_define_flag();
 
152
    parm.no_solpos->key = 's';
 
153
    parm.no_solpos->description = _("Do not use SOLPOS algorithm of NREL");
 
154
 
 
155
    if (G_parser(argc, argv))
 
156
        exit(EXIT_FAILURE);
 
157
 
 
158
    G_get_window(&window);
 
159
 
 
160
    /* require at least one output */
 
161
    elev_name = parm.elev->answer;
 
162
    azimuth_name = parm.azimuth->answer;
 
163
    sunhour_name = parm.sunhours->answer;
 
164
    if (!elev_name && !azimuth_name && !sunhour_name)
 
165
        G_fatal_error(_("No output requested, exiting."));
 
166
 
 
167
    year = atoi(parm.year->answer);
 
168
    if (parm.month->answer)
 
169
        month = atoi(parm.month->answer);
 
170
    else
 
171
        month = -1;
 
172
 
 
173
    day = atoi(parm.day->answer);
 
174
    hour = atoi(parm.hour->answer);
 
175
    minutes = atoi(parm.minutes->answer);
 
176
    seconds = atoi(parm.seconds->answer);
 
177
 
 
178
    lst_time = (parm.lst_time->answer != 0);
 
179
    use_solpos = (parm.no_solpos->answer == 0);
 
180
 
 
181
    /* init variables */
 
182
    ha = 180;
 
183
    ha_cos = 0;
 
184
    sd_cos = 0;
 
185
    sd_sin = 1;
 
186
    north_gc_cos = 0;
 
187
    north_gc_sin = 1;
 
188
    se_sin = 0;
 
189
 
 
190
    if (use_solpos && lst_time) {
 
191
        G_warning(_("NREL SOLPOS algorithm uses Greenwich standard time."));
 
192
        G_warning(_("Time will be interpreted as Greenwich standard time."));
 
193
 
 
194
        lst_time = 0;
 
195
    }
 
196
    if (!use_solpos) {
 
197
        if (lst_time)
 
198
            G_message(_("Time will be interpreted as local sidereal time."));
 
199
        else
 
200
            G_message(_("Time will be interpreted as Greenwich standard time."));
 
201
        
 
202
        if (sunhour_name)
 
203
            G_fatal_error(_("Sunshine hours require NREL SOLPOS."));
 
204
    }
 
205
 
 
206
    if ((G_projection() != PROJECTION_LL)) {
 
207
        if (window.proj == 0)
 
208
            G_fatal_error(_("Current projection is x,y (undefined)."));
 
209
 
 
210
        do_reproj = 1;
 
211
 
 
212
        /* read current projection info */
 
213
        if ((in_proj_info = G_get_projinfo()) == NULL)
 
214
            G_fatal_error(_("Cannot get projection info of current location"));
 
215
 
 
216
        if ((in_unit_info = G_get_projunits()) == NULL)
 
217
            G_fatal_error(_("Cannot get projection units of current location"));
 
218
 
 
219
        if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
 
220
            G_fatal_error(_("Cannot get projection key values of current location"));
 
221
 
 
222
        G_free_key_value(in_proj_info);
 
223
        G_free_key_value(in_unit_info);
 
224
 
 
225
        /*  output projection to lat/long w/ same ellipsoid as input */
 
226
        oproj.zone = 0;
 
227
        oproj.meters = 1.;
 
228
        sprintf(oproj.proj, "ll");
 
229
        if ((oproj.pj = pj_latlong_from_proj(iproj.pj)) == NULL)
 
230
            G_fatal_error(_("Unable to update lat/long projection parameters"));
 
231
    }
 
232
 
 
233
    /* always init pd */
 
234
    S_init(&pd);
 
235
    pd.function = S_GEOM;
 
236
    if (use_solpos) {
 
237
        pd.function = S_ZENETR;
 
238
        if (azimuth_name)
 
239
            pd.function = S_SOLAZM;
 
240
        if (sunhour_name)
 
241
            pd.function |= S_SRSS;
 
242
    }
 
243
    if (month == -1)
 
244
        doy = day;
 
245
    else
 
246
        doy = dom2doy2(year, month, day);
 
247
    
 
248
    set_solpos_time(&pd, year, 1, doy, hour, minutes, seconds);
 
249
    set_solpos_longitude(&pd, 0);
 
250
    pd.latitude = 0;
 
251
    S_solpos(&pd);
 
252
 
 
253
    if (lst_time) {
 
254
        /* hour angle */
 
255
        /***************************************************************
 
256
         * The hour angle of a point on the Earth's surface is the angle
 
257
         * through which the earth would turn to bring the meridian of
 
258
         * the point directly under the sun. This angular displacement
 
259
         * represents time (1 hour = 15 degrees).
 
260
         * The hour angle is negative in the morning, zero at 12:00,
 
261
         * and positive in the afternoon
 
262
         ***************************************************************/
 
263
 
 
264
        ha = 15.0 * (hour + minutes / 60.0 + seconds / 3600.0) - 180.;
 
265
        G_debug(1, "Solar hour angle, degrees: %.2f", ha);
 
266
        ha *= DEG2RAD;
 
267
        ha_cos = cos(ha);
 
268
        roundoff(&ha_cos);
 
269
    }
 
270
 
 
271
    if (!use_solpos) {
 
272
        /* sun declination */
 
273
        /***************************************************************
 
274
         * The declination of the sun is the angle between
 
275
         * the rays of the sun and the plane of the Earth's equator.
 
276
         ***************************************************************/
 
277
 
 
278
        s_gamma = (2 * M_PI * (doy - 1)) / 365;
 
279
        G_debug(1, "fractional year in radians: %.2f", s_gamma);
 
280
        /* sun declination for day of year with Fourier series representation
 
281
         * NOTE: based on 1950, override with solpos */
 
282
        s_declination = (0.006918 - 0.399912 * cos(s_gamma) + 0.070257 * sin(s_gamma) -
 
283
                         0.006758 * cos(2 * s_gamma) + 0.000907 * sin(2 * s_gamma) -
 
284
                         0.002697 * cos(3 * s_gamma) + 0.00148 * sin(3 * s_gamma));
 
285
 
 
286
        G_debug(1, "sun declination: %.5f", s_declination * RAD2DEG);
 
287
        G_debug(1, "sun declination (solpos): %.5f", pd.declin);
 
288
 
 
289
        if (lst_time) {
 
290
            north_ll = (window.north + window.south) / 2;
 
291
            east_ll = (window.east + window.west) / 2;
 
292
            if (do_reproj) {
 
293
                if (pj_do_proj(&east_ll, &north_ll, &iproj, &oproj) < 0)
 
294
                    G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
 
295
            }
 
296
            pd.timezone = east_ll / 15.;
 
297
            pd.time_updated = 1;
 
298
            set_solpos_longitude(&pd, east_ll);
 
299
            G_debug(1, "fake timezone: %.2f", pd.timezone);
 
300
            S_solpos(&pd);
 
301
            G_debug(1, "Solar hour angle (solpos), degrees: %.2f", pd.hrang);
 
302
        }
 
303
 
 
304
        /* always use solpos sun declination */
 
305
        s_declination = pd.declin * DEG2RAD;
 
306
        sd_sin = sin(s_declination);
 
307
        roundoff(&sd_sin);
 
308
        sd_cos = cos(s_declination);
 
309
        roundoff(&sd_cos);
 
310
 
 
311
        G_debug(1, "sun declination (solpos): %.5f", s_declination * RAD2DEG);
 
312
 
 
313
        if (0 && lst_time) {
 
314
            pd.timezone = 0;
 
315
            pd.time_updated = pd.longitude_updated = 1;
 
316
            S_solpos(&pd);
 
317
        }
 
318
    }
 
319
 
 
320
    if (elev_name) {
 
321
        if ((elev_fd = Rast_open_new(elev_name, FCELL_TYPE)) < 0)
 
322
            G_fatal_error(_("Unable to create raster map <%s>"), elev_name);
 
323
 
 
324
        elevbuf = Rast_allocate_f_buf();
 
325
    }
 
326
    else {
 
327
        elevbuf = NULL;
 
328
        elev_fd = -1;
 
329
    }
 
330
 
 
331
    if (azimuth_name) {
 
332
        if ((azimuth_fd = Rast_open_new(azimuth_name, FCELL_TYPE)) < 0)
 
333
            G_fatal_error(_("Unable to create raster map <%s>"), azimuth_name);
 
334
 
 
335
        azimuthbuf = Rast_allocate_f_buf();
 
336
    }
 
337
    else {
 
338
        azimuthbuf = NULL;
 
339
        azimuth_fd = -1;
 
340
    }
 
341
 
 
342
    if (sunhour_name) {
 
343
        if ((sunhour_fd = Rast_open_new(sunhour_name, FCELL_TYPE)) < 0)
 
344
            G_fatal_error(_("Unable to create raster map <%s>"), sunhour_name);
 
345
 
 
346
        sunhourbuf = Rast_allocate_f_buf();
 
347
    }
 
348
    else {
 
349
        sunhourbuf = NULL;
 
350
        sunhour_fd = -1;
 
351
    }
 
352
 
 
353
    if (elev_name && azimuth_name) {
 
354
        G_message(_("Calculating solar elevation and azimuth..."));
 
355
    }
 
356
    else if (elev_name) {
 
357
        G_message(_("Calculating solar elevation..."));
 
358
    }
 
359
    else if (azimuth_name) {
 
360
        G_message(_("Calculating solar azimuth..."));
 
361
    }
 
362
 
 
363
    nrows = Rast_window_rows();
 
364
    ncols = Rast_window_cols();
 
365
 
 
366
    ba2 = 6356752.3142 / 6378137.0;
 
367
    ba2 = ba2 * ba2;
 
368
 
 
369
    for (row = 0; row < nrows; row++) {
 
370
 
 
371
        G_percent(row, nrows, 2);
 
372
        
 
373
        /* get cell center northing */
 
374
        north = window.north - (row + 0.5) * window.ns_res;
 
375
        north_ll = north;
 
376
 
 
377
        for (col = 0; col < ncols; col++) {
 
378
            long int retval;
 
379
            s_elevation = s_azimuth = -1.;
 
380
 
 
381
            /* get cell center easting */
 
382
            east = window.west + (col + 0.5) * window.ew_res;
 
383
            east_ll = east;
 
384
 
 
385
            if (do_reproj) {
 
386
                north_ll = north;
 
387
 
 
388
                if (pj_do_proj(&east_ll, &north_ll, &iproj, &oproj) < 0)
 
389
                    G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
 
390
            }
 
391
 
 
392
            /* geocentric latitude */
 
393
            north_gc = atan(ba2 * tan(DEG2RAD * north_ll));
 
394
            north_gc_sin = sin(north_gc);
 
395
            roundoff(&north_gc_sin);
 
396
            north_gc_cos = cos(north_gc);
 
397
            roundoff(&north_gc_cos);
 
398
 
 
399
            if (!lst_time) {
 
400
                set_solpos_longitude(&pd, east_ll);
 
401
                pd.latitude = north_gc * RAD2DEG;
 
402
                retval = S_solpos(&pd);
 
403
                S_decode(retval, &pd);
 
404
                G_debug(3, "solpos hour angle: %.5f", pd.hrang);
 
405
            }
 
406
 
 
407
            /* solar elevation angle */
 
408
            if (!use_solpos) {
 
409
                if (!lst_time) {
 
410
                    ha = pd.hrang;
 
411
                    ha_cos = cos(ha * DEG2RAD);
 
412
                    roundoff(&ha_cos);
 
413
                }
 
414
                se_sin = ha_cos * sd_cos * north_gc_cos + sd_sin * north_gc_sin;
 
415
                roundoff(&se_sin);
 
416
                s_elevation = RAD2DEG * asin(se_sin);
 
417
            }
 
418
            else /* use_solpos && !lst_time */
 
419
                s_elevation = pd.elevetr;
 
420
 
 
421
            if (elev_name)
 
422
                elevbuf[col] = s_elevation;
 
423
 
 
424
            if (azimuth_name) {
 
425
                /* solar azimuth angle */
 
426
                if (!use_solpos) {
 
427
                    sa_cos = (se_sin * north_gc_sin - sd_sin) /
 
428
                             (cos(DEG2RAD * s_elevation) * north_gc_cos);
 
429
                    roundoff(&sa_cos);
 
430
                    s_azimuth = RAD2DEG * acos(sa_cos);
 
431
                    
 
432
                    /* morning */
 
433
                    s_azimuth = 180. - RAD2DEG * acos(sa_cos);
 
434
                    if (ha > 0)   /* afternoon */
 
435
                        s_azimuth = 360.0 - s_azimuth;
 
436
                }
 
437
                else
 
438
                    s_azimuth = pd.azim;
 
439
 
 
440
                azimuthbuf[col] = s_azimuth;
 
441
            }
 
442
            
 
443
            if (sunhour_name) {
 
444
                sunhourbuf[col] = (pd.ssetr - pd.sretr) / 60.;
 
445
                if (sunhourbuf[col] > 24.)
 
446
                    sunhourbuf[col] = 24.;
 
447
                if (sunhourbuf[col] < 0.)
 
448
                    sunhourbuf[col] = 0.;
 
449
            }
 
450
 
 
451
        }
 
452
        if (elev_name)
 
453
            Rast_put_f_row(elev_fd, elevbuf);
 
454
        if (azimuth_name)
 
455
            Rast_put_f_row(azimuth_fd, azimuthbuf);
 
456
        if (sunhour_name)
 
457
            Rast_put_f_row(sunhour_fd, sunhourbuf);
 
458
    }
 
459
    G_percent(1, 1, 2);
 
460
 
 
461
    if (elev_name) {
 
462
        Rast_close(elev_fd);
 
463
        /* writing history file */
 
464
        Rast_short_history(elev_name, "raster", &hist);
 
465
        Rast_command_history(&hist);
 
466
        Rast_write_history(elev_name, &hist);
 
467
    }
 
468
    if (azimuth_name) {
 
469
        Rast_close(azimuth_fd);
 
470
        /* writing history file */
 
471
        Rast_short_history(azimuth_name, "raster", &hist);
 
472
        Rast_command_history(&hist);
 
473
        Rast_write_history(azimuth_name, &hist);
 
474
    }
 
475
    if (sunhour_name) {
 
476
        Rast_close(sunhour_fd);
 
477
        /* writing history file */
 
478
        Rast_short_history(sunhour_name, "raster", &hist);
 
479
        Rast_command_history(&hist);
 
480
        Rast_write_history(sunhour_name, &hist);
 
481
    }
 
482
 
 
483
    G_done_msg(" ");
 
484
 
 
485
    exit(EXIT_SUCCESS);
 
486
}
 
487
 
 
488
void set_solpos_time(struct posdata *pdat, int year, int month, int day,
 
489
                    int hour, int minute, int second)
 
490
{
 
491
    pdat->year = year; 
 
492
    pdat->month = month; 
 
493
    pdat->day = day; 
 
494
    pdat->daynum = day; 
 
495
    pdat->hour = hour; 
 
496
    pdat->minute = minute; 
 
497
    pdat->second = second;
 
498
    pdat->timezone = 0;
 
499
 
 
500
    pdat->time_updated = 1; 
 
501
    pdat->longitude_updated = 1;
 
502
}
 
503
 
 
504
void set_solpos_longitude(struct posdata *pdat, double longitude)
 
505
{
 
506
    pdat->longitude = longitude;
 
507
 
 
508
    pdat->longitude_updated = 1;
 
509
}
 
510
 
 
511
int roundoff(double *x)
 
512
{
 
513
    /* watch out for the roundoff errors */
 
514
    if (fabs(*x) > 1.0) {
 
515
        if (*x > 0.0)
 
516
            *x = 1.0;
 
517
        else
 
518
            *x = -1.0;
 
519
 
 
520
        return 1;
 
521
    }
 
522
 
 
523
    return 0;
 
524
}