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