GRASS-Wiki - User contributions [en]

V.generalize tutorial

2026-03-18T05:02:48Z

Hamish: /* Simplification */ add some explaining

== MODULE ==

* GRASS 7.x {{cmd|v.generalize}}

== TUTORIAL ==

== Introduction ==

This tutorial presents and explains the functionality of GRASS vector module ''v.generalize''.
This module implements generalization operators for GRASS vector maps. The topics covered in this tutorial are: simplification, smoothing, network generalization and displacement. For the basic introduction to these operators check the {{cmd|v.generalize}} man page, or for the exhausting introduction check McMaster and Shea (TODO: add reference).

It is recommended to read the official man page before reading this document, or read
these two documents at the same time since this tutorial does not contain the
(detail) descriptions of input parameters whereas the man page has very few examples
and no pictures at all. (Even the author of module and both of the documents must have consult the man page several times....)
This document is also meant to be a "report" showing the work done during the Google Summer of Code 2007.

Most of the examples in this text work with Spearfish dataset, which can
be downloaded [[Sample_datasets|here]]. Also, this
tutorial assumes that the user is already running GRASS session with the Spearfish
location and a monitor opened. Also, if you click on any picture here, it will be shown in its full size.

All the algorithms presented in this document (try to) preserve the topology
of the input maps. This means, for example, that the smoothing and simplification methods
never remove the first and last points of lines and that displacement preserves the junctions.

== Simplification ==

''v.generalize'' implements many simplification algorithms. The most widely used
is Douglas-Peucker algorithm (TODO: reference). We can apply this algorithm to
any vector map in the following way:

<pre>
v.generalize input=roads output=roads_douglas method=douglas threshold=50
</pre>

The explanation of the line above is following: Run v.generalize, apply Douglas-Peucker
algorithm with threshold equals to 50 to the map roads and store the output in vector
map output_douglas. The last output line produced by this module:

<pre>
Number of vertices was reduced from 5468 to 2107[38%]
</pre>

means that the input file (roads) has 5468 vertices in total and the new map (roads_douglas)
has only 2107 vertices which is only 38% of original. On the other hand, if we run the
commands:

<pre>
d.vect roads
d.vect roads_douglas col=blue
</pre>

We can see that there are no significant differences between the input and output maps. "Only the details were removed".

[[Image:v.generalize.pic1.png|center|thumb|400px|caption text here]]

The amount of the details removed can be specified by parameter: threshold. It is the
case that the output map has fewer vertices and details for greater values of threshold.
For example, if we run

<pre>
v.generalize input=roads output=roads_douglas method=douglas threshold=200 --overwrite
</pre>

we obtain a map with only 1726 vertices. A disadvantage of the command above is that
it never removes the lines. <strike>If we also want to remove the small lines, we need to run
the command above with the -r flag. If the -r flag is presented, lines shorter than
threshold and areas with areas less than threshold are removed:

<pre>
v.generalize -r input=roads output=roads_douglas method=douglas threshold=200 --overwrite
</pre>
</strike>
Fixme: -r flag and method=remove_small have been removed. Try v.clean for removing small areas or
v.edit map=mapname type=line tool=delete query=length thresh=0,0,-0.5
<strike>
In this case, roads_douglas has only 850 vertices and it contains 387 lines whereas
the original map (roads) contains 825 lines. In this case, the output map has very few details,
but the basic shapes and topology are preserved:

[[Image:v.generalize.pic2.png|center|thumb|400px|caption text here]]
</strike>
<strike>It is also possible remove small lines/areas only (without any simplification). This is achieved
by method=remove_small:

<pre>
v.generalize input=roads output=roads_remove_small method=remove_small threshold=200
</pre>

produces following map with 399 lines. (Note removed lines in the upper left corner)

[[Image:v.generalize.pic3.png|center|thumb|400px|caption text here]]
</strike>

Douglas-Peucker Algorithm has very reasonable results, but it is very hard to find (i.e. guess) the
correct value to use for ''threshold=''. Moreover, it is also impossible to simplify each line to (for example) 40%.
Exactly for such cases, ''v.generalize'' provides method=douglas_reduction. This algorithm
is a modification of Douglas-Peucker Algorithm which takes another parameter ''reduction''
which denotes (approximate, +/- a few percent) simplification of lines. A value of ''reduction=10'' will keep (approximately) 10% of the original vertices and drop the remaining 90%.

For example,

<pre>
v.generalize input=roads output=roads_douglas_reduction method=douglas_reduction \
threshold=0 reduction=50 --overwrite
d.erase
d.vect roads_douglas_reduction
</pre>

produces following map with 3018 vertices (55%). (Note that there are almost no differencies between the original and the new map)

[[Image:v.generalize.pic4.png|center|thumb|400px|caption text here]]

Also observe that the following commands are equivalent

<pre>
v.generalize input=in output=out method=douglas threshold=eps
v.generalize input=in output=out method=douglas_reduction threshold=eps reduction=100
</pre>

Another algorithm implemented in this modules is "Vertex Reduction". This algorithm removes
the consecutive poins (on the same line) which are closer to each other than threshold. For example,

<pre>
v.generalize input=in output=out method=reduction threshold=0
</pre>

removes duplicate points. More precisely, if two consecutive points have the same
coordinates then the second point is removed and the first is preserved. The last two algorithm
implemented by this module are "Lang" and "Reumann-Witkam" algorithm. For more information about these two algorithms, please see the {{cmd|v.generalize}} man page.

== Results ==

The following four pictures show the results obtained by Reumann, Douglas, Lang
and Vertex Reduction algorithm resp. The algorithms were run with threshold set to 50
and Lang algorithm with look_ahead=7.

<gallery Caption="Simplification algorithms" widths="300px" heights="225px" perrow="2">
Image:v.generalize.reumann.png|Reumann-Witkam algorithm result containing 2522 [46%] points

Image:v.generalize.douglas.png|Douglas algorithm result containing 2107 [38%] points

Image:v.generalize.lang.png|Lang algorithm result containing 2160 [39%] points

Image:v.generalize.reduction.png|Vertex Reduction algorithm result containing 4296 [78%] points
</gallery>

The map produced by
* Reumann-Witkam algorithm contains 2522 [46%] points,
* Douglas: 2107 [38%] points,
* Lang: 2160 [39%] and
* Vertex Reduction: 4296 [78%].

== Smoothing ==

''v.generalize'' also supports many smoothing algorithm. For basic descriptions, please consult the {{cmd|v.generalize}} man page.

Probably, the best results are produced by "Chaiken", "Hermite" and "Snakes" algorithms/methods.
However, the remaining algorithms may also produce very reasonable results. Although
the Chaiken and Hermite methods may produce the maps with a lot of new points, the methods
presented above (simplification) provide a good tool for tackling this problem.

If we run the following command
<pre>
v.generalize input=roads output=roads_chaiken method=chaiken threshold=1
</pre>

we get a new map with 33364[640%] vertices.

[[Image:v.generalize.pic5.png|center|thumb|400px|caption text here]]

This map looks almost exactly the same as the original map at the current level of detail as the
picture below shows. This pictures was produced by the command:

<pre>
d.erase
d.vect roads
d.vect roads_chaiken col=blue
</pre>

[[Image:v.generalize.pic6.png|center|thumb|400px|caption text here]]

However, if we zoom to a small region, we can see that the new map consists of smooth(er) lines
which (very reasonable) approximate the original ones.

[[Image:v.generalize.pic7.png|center|thumb|400px|caption text here]]

If we apply "Hermite" method instead, we will obtain a map with 14640[267%] vertices.

<pre>
v.generalize input=roads output=roads_hermite method=hermite threshold=1
</pre>

[[Image:v.generalize.pic8.png|center|thumb|400px|caption text here]]

Note, that a difference between "Chaiken" and "Hermite" is that the lines produced
by "Chaiken" "inscribe" the original lines whereas the "Hermite" lines "circumscribe" the
original lines as can be seen in the picture below. (Black line is original line,
green line is "Chaiken" and blue is "Hermite")

[[Image:v.generalize.pic9.png|center|thumb|400px|caption text here]]

The algorithms mentioned above are suitable for smooth approximation of given lines.
On the other hand, if the aim of smoothing is to "straighten" the lines then the better
results are achieved by the other methods. For example,

<pre>
v.generalize input=roads output=roads_sa method=sliding_averaging look_ahead=7 slide=1
</pre>



At first sight, we can see that roads_sa contains smooth and straight lines which
preserve the original shape of the lines. This difference is obvious if we zoom to a small
region of a map (see below. Again, original line is black, new line is blue)

[[Image:v.generalize.pic11.png|center|thumb|400px|caption text here]]

If the lines are "too straight" then we can set "slide" to a smaller value to obtain the lines
which better preserve the original shape. In the picture below, original line is black,
line produced by "slide=1" is blue and "slide=0.3" is green.

[[Image:v.generalize.pic12.png|center|thumb|400px|caption text here]]

Very similar results can be obtained by Distance Weighting Algorithm (method=distance_weighting).
This is not very surprising since these algorithms are almost the same. For example,
the image below shows the outputs of "Distance Weighting Algorithm". The image was generated
by the following sequence of commands

<pre>
v.generalize input=roads output=roads_dw1 method=distance_weighting look_ahead=7 slide=1
v.generalize input=roads output=roads_dw2 method=distance_weighting look_ahead=7 slide=0.3
d.erase
d.vect roads
d.vect roads_dw1 col=red
d.vect roads_dw2 col=blue
</pre>

[[Image:v.generalize.pic13.png|center|thumb|400px|method=distance_weighting]]

Also, very good resutls can be obtained by the "Snakes" algorithm. On the other hand,
it is the (asymptotically) slowest
smoothing algorithm implemented in this however. Behaviour of this algorithm is controlled by
"alpha" and "beta" parameter. Reasonable range of values for these two parameters is [0..5]
where higher values correspond to the straighter lines. Module outputs the input map if
alpha=beta=0. And this command

<pre>
v.generalize input=roads output=roads_snakes method=snakes alpha=1 beta=1
</pre>

produces a map containing following region (original line is black)
[[Image:v.generalize.pic14.png|center|thumb|400px|caption text here]]

Last smoothing algorithm implemented in this module is "Boyle's Forward-Looking Algorithm" which
is another "straightening" algorithm.

<pre>
v.generalize input=roads output=roads_boyle method=boyle look_ahead=5
</pre>

produces a map containing following region (original line is black)
[[Image:v.generalize.pic15.png|center|thumb|400px|caption text here]]

=== Area smoothing example ===

# spearfish
g.region rast=geology
r.reclass in=geology out=geology.claysand << EOF
8 = 8 claysand
EOF
r.to.vect in=geology.claysand out=geology_claysand feature=area
v.generalize in=geology_claysand out=geology_claysand_smooth method=snakes

== Displacement ==

If we render entire Spearfish location, we can see in the upper half of the map two interstates which
overlap. This is not logically correct (I hope, they do not evarlap in real) and it is also
considered as an (presentation) error. For solving such problems, ''v.generalize'' provides
"dislplacement" method. As the name suggests, this method displaces linear features which are
close to each other so that they do not overlap/collide. Method implemented in this modules (based on Snakes)
has very good results but not very good perfomance. Therefore the calculations may take few(several) minutes.
For this reason, displacement is applied to the simplified lines in this document.

<pre>
v.generalize input=roads output=roads_dr method=douglas_reduction threshold=0 reduction=50
v.generalize input=roads_dr output=roads_dr_disp method=displacement alpha=0.01 beta=0.01 threshold=100 iterations=35
</pre>

First command produces simplified lines and then the second command applied displacement operator
to the simplified line. Parameters alpha and beta specifies the rigidity of the lines. This means
that displacement is bigger for small values of alpha and beta. Also, the displacement is
not very significant for higher(>=1.0) values of alpha, beta. Threshold parameter
denotes the critical distance. Only the points (and their neighbours) which are closer
than threshold apart are displaced by ''v.generalize''. Module tries to move these
points such that they are at least threshold apart. However, the displaced points are never
threshold (or more) apart for positive values of alpha and beta. Displacement as implemented
in ''v.generalize'' is an iterative process. Parameter "iterarions" specifies the number
of iterations the collisions between the points are resolved. In general, the quality of displacement
increases with the number of iterations. However, quality converges quite rapidly and
for all maps I tried, the sufficient value of iterations was between 20 and 50.
Two command above produce the picture below. Note that it is now possible to distinguish
between two "interastate lines" and also observe the free space between interstate and the lines
directly below it.

[[Image:v.generalize.pic16.png|center|thumb|400px|caption text here]]

== Network Generalization ==

Network generalization is suitable for selecting "the most important" subnetwork
of a network. For example, to select highways, interstates from a road network. Examples
in this section work with new GRASS default dataset, which can be downloaded [[Sample_datasets|here]].

If we render map "streets_wake" we really cannot see the streets, but the only
thing we can see is a big black rhombus. We will try to improve this. Firstly, network
generalization requites quite a lot of time and memory. Therefore, we begin with
simplification of "streets_wake":

<pre>
v.generalize input=streets_wake output=streets_rs method=remove_small threshold=50
</pre>

Then we can begin with network generalization. If we execute the folllowing command:

<pre>
v.generalize input=streets_rs output=streets_rs_network method=network betweeness_thresh=50
</pre>

we obtain the following map containing "only" 14128. Original map has 49746

[[Image:v.generalize.net50.png|center|thumb|400px|caption text here]]

If this is still not enough, we can increase the value of betweeness_thresh to, for example 200.
For this value, ''v.generalize'' produces following map with 11537 lines.

[[Image:v.generalize.net200.png|center|thumb|400px|caption text here]]

It is also possible to change the values of "closeness_thresh" and "degree_thresh". Parameter
"closeness_thresh" is suitable for selecting the "centre(s)" of a network. This parameter
is always between 0 and 1. And "reasonable values" of this parameter are smaller for
bigger network.

== Gereneral Parameters ==

''v.generalize'' has some parameters and flags which affect the general behaviour of module.

The simplest one is -c flag. "C" stands for copy and if this flag is on then the attributes
are copied from the old map to the new map. Note that the attributes of removed features are
dropped.

Default behaviour of this module is that the selected algorithm/method is applied to the all
lines/areas. It is possible to apply the most of the algorithms only to the selected features.
This is achieved by "type", "layer", "cats" and "where" parameters. This works for all algorithms
except "Network Generalization" which is always applied to the all features.
For example, the following command applies "Douglas Reduction" algorithm to interstates and highways (cat<3)
and leaves the other lines unaltered. It also copies the attributes

<pre>
v.generalize -c input=roads output=roads_douglas_reduction2 method=douglas_reduction threshold=0 reduction=50 type=line layer=1 where="cat<3"
</pre>

And the following command removes the small areas

<pre>
v.generalize input=soils output=soils_remove_small method=remove_small threshold=200 type=area
</pre>

Similarly, the following command displaces only the interstates (cats=1) and the lines with a different
category number are not taken into the account.

<pre>
v.generalize input=roads output=roads_displacement2 method=displacement \
threshold=75 alpha=0.01 beta=0.01 iterations=20 layer=1 cats=1
</pre>

[[Image:v.generalize.pic10.png|center|thumb|400px|caption text here]]

We end up with a complex example of a generalization of "roads" in Spearfish location.

<pre>
#straighten the lines
v.generalize input=roads output=step1 method=snakes alpha=1 beta=1
#simplification
v.generalize input=step1 output=step2 method=douglas_reduction threshold=0 reduction=55
#displacement
v.generalize input=step2 output=step3 method=displacement alpha=0.01 beta=0.01 threshold=100 iterations=50
#remove small areas
v.generalize input=step3 output=step4 method=remove_small threshold=75
#network generalization
v.generalize input=step4 output=step5 method=network betweeness_thresh=5 closeness_thresh=0.0425
#smoothing
v.generalize input=step5 output=step6 method=chaiken threshold=1
#simplification
v.generalize input=step6 output=step7 method=douglas threshold=1
#remove temporary maps
g.remove vect=step1,step2,step3,step4,step5,step6
</pre>

The result "step7" has 655 lines and 3545 vertices and the commands above have the following effect:

[[Image:v.generalize.anim.gif|center|thumb|400px|caption text here (animation)]]

== References ==

TODO

== See also ==

Use cases:
* http://www.custom-scenery.org/Line-data-simpl.275.0.html

== AUTHORS ==

* Daniel Bundala, Google Summer of Code 2007, Student
* Wolf Bergenheim, Mentor

[[Category: Documentation]]
[[Category: Tutorial]]
[[Category: Vector]]

Movies

2025-11-04T03:18:36Z

Hamish: /* Encode as animated GIF */

==== Q: How to construct high-quality MPEG-4 movies (animations) from a series of still frames? ====

''Maris wrote'':

Encode all .png files in directory to out.avi with 15 frames per second:
<pre>
mencoder "mf:///path/to/files/*.png" -mf fps=15 -o out.avi # put -ovc here, see next lines
</pre>
<pre>
## for DivX - libavcodec MPEG 4 (DivX5), use:
-ovc lavc -lavcopts vcodec=mpeg4:vhq:vbitrate=1800 -ffourcc MP4S
## for XviD, use:
-ovc xvid -xvidencopts bitrate=1024 -xvidencopts pass=2 -xvidencopts cartoon
## for DivX4, use:
-ovc divx4 -divx4opts q=5
## for x264 use:
-ovc x264 -x264encopts pass=3 -x264encopts preset=veryslow -x264encopts tune=animation
</pre>

Uses mplayer's encoder. Choose one of encoding formats and append to end
of "mencoder" line. Unfortunately only way to know quality settings is by
encoding, watching result and re-encoding with different params.

Notice - by default encoding with lavc will set video fourcc to FMP4, that can be decoded with ffdshow. Option -ffourcc MP4S will change it to MPS4 ("official" MS fourcc for mpeg4) and video will be playable on Windows by standard MS mpeg4 decoder. It may lead to side effects, if MS mpeg4 decoder is buggy. More info: http://www.itdp.de/mplayer-users/2005-03/msg00069.html

Other guides:
* http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-enc-images.html
* http://electron.mit.edu/~gsteele/ffmpeg/
* http://www.stillhq.com/extracted/howto-jpeg2mpeg/output.html

--------

==== Q: How to construct high-quality WebM movies? ====

A: {{wikipedia|WebM}} is the container for the {{wikipedia|VP8}} and {{wikipedia|VP9}} codecs. VP8 approaches H.264 in quality but differs from H.264 in that it is both patent + license free to use. Likewise, VP9 is an alternative to HEVC/H.265. The audio channel, if used, is encoded using the Ogg {{wikipedia|Vorbis}} codec for VP8 or Opus codec for VP9. Google has been promoting WebM as the standard fallback video format for HTML5 and is moving YouTube content to it. Firefox. Opera, and Google Chrome web browsers support WebM playback out of the box without need for extra plugins.

; Encoding:

First convert PNG or JPEG images to PPM format using the {{wikipedia|NetPBM}} tools.
: ''This step is pretty fast, I wouldn't worry about trying to parallelize it.''

<syntaxhighlight lang="bash">
for file in *.png ; do
pngtopnm $file > `basename $file .png`.ppm
done
</syntaxhighlight>
::or
<syntaxhighlight lang="bash">
for file in *.jpg ; do
jpegtopnm $file > `basename $file .jpg`.ppm
done
</syntaxhighlight>

Next create a Y4M stream containing a concatenation of the raw PPM frames:

<syntaxhighlight lang="bash">
cat *.ppm | ppmtoy4m -F 2:1 > all_frames.y4m
</syntaxhighlight>

* ''The "-F 2:1" option above sets the frame rate to 2 frames per second.''
* ''The <tt>ppmtoy4m</tt> program comes from the mjpegtools package.''
* ''For a very high number of frames you may wish to skip this step and create a named FIFO or stdin pipe directly into <tt>vpxenc</tt>, as the raw Y4M file can become extremely large.''

Finally, run the <tt>vpxenc</tt> encoder to create a high quality WebM video:
<syntaxhighlight lang="bash">
vpxenc all_frames.y4m -o my_movie.webm \
--passes=2 --threads=4 \
--best --target-bitrate=2000 --end-usage=vbr \
--auto-alt-ref=1 --verbose \
--minsection-pct=5 --maxsection-pct=800 \
--lag-in-frames=16 --kf-min-dist=0 --kf-max-dist=360 \
--static-thresh=0 --drop-frame=0 \
--min-q=0 --max-q=60
</syntaxhighlight>

* ''Set "--threads=" in the above command to the number of cores on your CPU. (or one less than the number of cores if you have many)''
* ''If you are encoding to VP8 (and not VP9) consider adding <tt>--token-parts=2</tt>''
* ''If <tt>vpxenc</tt> complains about "Invalid image format" either add "-S 420jpeg" to the <tt>ppmtoy4m</tt> options or add "--profile=1" to the <tt>vpxenc</tt> options.
* ''The <tt>vpxenc</tt> program comes from the vpx-tools package.''
* [http://www.webmproject.org/docs/encoder-parameters/#2-pass-best-quality-vbr-encoding Encoder parameters help page]

--------

==== Q: How to generate an animated GIF, FLI, or MNG? ====

A: for animations of 300 frames or less, animated GIF, [http://woodshole.er.usgs.gov/operations/modeling/flc.html FLI], or MNG formats are smaller files and better quality (i.e., frames stored as PNG not JPEG). Create with [http://www.lcdf.org/gifsicle/ gifsicle], [http://vento.pi.tu-berlin.de/fli.html ppm2fli], or for MNG use [http://www.imagemagick.org/script/convert.php ImageMagick's 'convert']. With more than 300 frames the players usually have memory issues. Animated GIFs are playable in any web browser of course and also in OpenOffice.

# with convert, many GIFs into one animated GIF:
convert -verbose -delay 20 -loop 0 snapshot*.png animation.gif

Notes: '-delay 20' means 20 hundreds of a second delay between each frame. The '-loop 0' flag lets it loop indefinitely. To loop only 3 times use '-loop 3' etc.

Sample script to export PNGs in GRASS (generates tmean_001.png .. tmean_708.png):

<syntaxhighlight lang=bash>
export i=0
for year in `seq 1950 2008` ; do
for m in `seq 1 12` ; do
i=`expr $i + 1`
# awk trick to generate wildcard-correct file names:
i=`echo $i | awk '{printf "%03d\n", $1}'`
r.out.png tmean_europe.$year.$m.avg out=tmean_$i.png
done
done
convert -verbose -delay 20 -loop 0 tmean_*.png tmean_animation.gif
</syntaxhighlight>

For presentations using a web browser, you can center the image on a blank white page, turn off any toolbars and go into full screen mode (F11 for Firefox or Chrome). Then Alt-Tab your way to the animation at the appropriate time and nobody knows you are using a web browser.

Example HTML for centering image:

<syntaxhighlight lang="html4strict">
<HTML>
<HEAD></HEAD>
<BODY>
 
<CENTER>
<IMG SRC="animation.gif">
</CENTER>
</BODY>
</HTML>
</syntaxhighlight>

--------

==== Q: How to create dynamic surface movies in NVIZ? ====

A: See slides from the FOSS4G 2006 workshop:
: http://skagit.meas.ncsu.edu/~helena/grasswork/foss4g/FOSS4G06WKSVisual4anim.odp

See also the example at the bottom of the [http://grass.osgeo.org/grass64/manuals/html64_user/nviz/nviz_panel_kanim.html NVIZ keyframe animator panel help page].

--------

==== Q: How to animate a time series of raster maps? ====

A: Use
* {{cmd|d.slide.show}} to view in a GRASS xmon (see {{cmd|d.mon}})
* {{cmd|r.out.mpeg}} to save to a MPEG-1 movie.
* GRASS GIS 6: {{cmd|xganim|version=64}} to view on screen in X-windows (with live start/stop/speed controls)
* GRASS GIS 7: {{cmd|g.gui.animation}} to view on screen in wxGUI (with live start/stop/speed controls)
: See also the wiki page for the [[wxGUI Animation Tool]]:
{{YouTube|Pi1ESk9iPwk|YouTube video introduction to the WxGUI Animation Tool}}

--------

==== Q: How to animate a series of maps, including decorations? ====

A: Use the xmon drivers ({{cmd|d.mon}})
to render the displays to a window via a shell script loop, then switch from the ''x0'' to ''PNG'' driver to save as a series of PNG or PPM images. Use a method answered in a FAQ above to combine those into a single animation.

--------

==== Q: How to create a screencast (for video tutorials etc)? ====

; GNU/Linux:
* Easist way: [http://recordmydesktop.sourceforge.net/ recordmydesktop] (user interfaces are 'qt-recordmydesktop' and 'gtk-recordmydesktop')
:* when uploading videos to Youtube created by <tt>recordmydesktop</tt> you need to encode the file eg. using <tt>mencoder</tt>
:- No sound, just video screencast, as AVI file:
::<pre>mencoder out.ogv -nosound -ovc lavc -o out.avi</pre>
:- With audio track, as AVI file:
::<pre>mencoder out.ogv -ovc xvid -oac mp3lame -xvidencopts pass=1 -o out.avi</pre>
:- With audio track, as WebM file:
::<pre>ffmpeg -i out.ogv -b:a 24k out.webm</pre>
:- With audio track, as MP4 file:
::<pre>ffmpeg -i out.ogv -ar 44100 -ab 96k -coder ac -strict experimental -vcodec libx264 out.mp4</pre>
* [https://wiki.ubuntu.com/MeetingLogs/openweekfeisty/screencast Ubuntu discussion over screencast]
* [https://wiki.ubuntu.com/ScreencastTeam/RecordingScreencasts Ubuntu screencast how to]
* [http://www.unixuser.org/~euske/vnc2swf/ vnc2swf] to record a flash movie ([http://www-pool.math.tu-berlin.de/~soeren/grass/modules/screenshots/vnc2swf_usage.html short vnc2swf usage tutorial])
* [http://www.maartenbaert.be/simplescreenrecorder/ SimpleScreenRecorder]

; MacOSX:
* [http://crschmidt.net/blog/archives/203/openlayers-screencast/ Chris Schmidt experience on MacOSX]

; MS Windows:
* [http://camstudio.org/ Camstudio]

-----

==== Examples ====
''Using the [http://grass.osgeo.org/download/data6.php Spearfish] sample dataset''

===== Animate on screen =====
''Using a shell script''

<syntaxhighlight lang=bash>
GRASS_WIDTH=640 GRASS_HEIGHT=500 \
d.mon start=x0

g.region rast=elevation.dem n=4927830 s=4912980

for NUM in `seq 1 0.5 10` ; do
d.rast elevation.dem --quiet
d.vect roads
d.vect bugsites icon=basic/circle col=black fcol=green size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph

echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
done
</syntaxhighlight>

===== Animate to files =====
''Using [http://grass.osgeo.org/grass64/manuals/html64_user/pngdriver.html GRASS's PNG driver]''

<syntaxhighlight lang=bash>
# set up base map
g.region rast=elevation.dem n=4927830 s=4912980
export GRASS_WIDTH=640
export GRASS_HEIGHT=500
export GRASS_PNGFILE=base.ppm
export GRASS_PNG_READ=FALSE
d.mon start=PNG
d.rast elevation.dem
d.vect roads
echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
d.mon stop=PNG

# loop to draw variable parts as new frames
i=0
export GRASS_PNG_READ=TRUE
for NUM in `seq 1 0.5 10` ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
echo "Processing frame $FRAMENUM ..."
export GRASS_PNGFILE="ganim_${FRAMENUM}.ppm"
cp base.ppm "$GRASS_PNGFILE"
d.mon start=PNG --quiet
d.vect bugsites icon=basic/circle col=black fcol=green \
size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph
d.mon stop=PNG --quiet
done

# repeat the last frame a few times
for EXTRA in 1 2 3 ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
cp base.ppm "ganim_${FRAMENUM}.ppm"
done
\rm base.ppm
</syntaxhighlight>

===== Encode as animated GIF =====

* ''Using the [http://www.lcdf.org/gifsicle/ gifsicle] encoder''
* ''View in a web browser or with ImageMagick's [http://www.imagemagick.org/Usage/basics/#animate animate]''
* see [[Movies#Q:_How_to_generate_an_animated_GIF.2C_FLI.2C_or_MNG.3F]] above

'''A) Using gifsicle:'''

Advantage: the delay is coded in the animated GIF file and rendered properly e.g. in a Web browser.

<syntaxhighlight lang=bash>
# Imagemagick's `convert` program may do a better job at
# 256 color quantization than NetPBM's `ppmquant`.
for IMG in ganim*.ppm ; do # convert each frame from PPM to GIF
#ppmquant 256 $IMG | ppmtogif > "`basename $IMG .ppm`.gif"
convert "$IMG" "`basename $IMG .ppm`.gif"
done

gifsicle -O2 --delay 20 --no-loopcount --colors 256 ganim_*.gif > ganim.gif

\rm ganim_0*.gif
</syntaxhighlight>

'''B) Using ImageMagick's convert and [http://www.imagemagick.org/Usage/basics/#animate animate]:'''
<syntaxhighlight lang=bash>
# Generate "film" (create animated GIF, 750 millisec between frames):
convert -delay 75 lst_monthly_avg_film.00*.png lst_monthly_avg_2002_2012_zoom_film.gif

# play animation (delay 500 millisec between frames, pause for 5 seconds at the end)
animate -delay 50 -pause 5 lst_monthly_avg_2002_2012_zoom_film.gif
</syntaxhighlight>

===== Encode as animated PNG (MNG) =====
* ''Using the [http://www.imagemagick.org ImageMagick] "convert" encoder''
* ''View with ImageMagick's "display"''

<pre>
convert -delay 20 ganim_*.ppm ganim.mng
</pre>

===== Encode as FLI =====
* ''Using the [http://vento.pi.tu-berlin.de/fli.html ppm2fli] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or [http://xanim.polter.net xanim]''

<syntaxhighlight lang=bash>
ls ganim_*.ppm > frames.txt
ppm2fli -g"${GRASS_WIDTH}x${GRASS_HEIGHT}" -O -s 15 frames.txt ganim.fli
\rm frames.txt
</syntaxhighlight>

===== Encode as Flash =====

* ''Using the [http://www.swftools.org/ SWFTools] <tt>png2swf</tt> encoder''
* ''View with [http://www.mplayerhq.hu mplayer]''

png2swf -o outfile.swf *.png -r 1

===== Encode as MPEG-4 Xvid =====
* ''Using the [http://www.mplayerhq.hu mencoder] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or similar''

<syntaxhighlight lang=bash>
# doesn't like GRASS's PPMs, convert to PNG
for IMG in ganim*.ppm ; do
pnmtopng $IMG > "`basename $IMG .ppm`.png"
done
mencoder "mf://ganim_*.png" -mf "type=png:fps=5" -o ganim.avi \
-ovc "xvid" -xvidencopts "bitrate=1024"
\rm ganim_0*.png
</syntaxhighlight>

==== Publish your cool videos in this Wiki ====
We have installed a Mediawiki widget extension and created a special template for including Youtube videos ([[Template:YouTube]]). See also [[WxGUI_Modeler#Video_tutorials|wxGUI Modeler]] page for the usage.

The overall goal is to have many as possible '''[http://www.youtube.com/results?search_query=grass+gis GRASS GIS videos on YouTube]'''

[[Category:FAQ]]
[[Category:Visualization]]

Movies

2025-11-04T03:18:16Z

Hamish: source -> syntaxhighlight

==== Q: How to construct high-quality MPEG-4 movies (animations) from a series of still frames? ====

''Maris wrote'':

Encode all .png files in directory to out.avi with 15 frames per second:
<pre>
mencoder "mf:///path/to/files/*.png" -mf fps=15 -o out.avi # put -ovc here, see next lines
</pre>
<pre>
## for DivX - libavcodec MPEG 4 (DivX5), use:
-ovc lavc -lavcopts vcodec=mpeg4:vhq:vbitrate=1800 -ffourcc MP4S
## for XviD, use:
-ovc xvid -xvidencopts bitrate=1024 -xvidencopts pass=2 -xvidencopts cartoon
## for DivX4, use:
-ovc divx4 -divx4opts q=5
## for x264 use:
-ovc x264 -x264encopts pass=3 -x264encopts preset=veryslow -x264encopts tune=animation
</pre>

Uses mplayer's encoder. Choose one of encoding formats and append to end
of "mencoder" line. Unfortunately only way to know quality settings is by
encoding, watching result and re-encoding with different params.

Notice - by default encoding with lavc will set video fourcc to FMP4, that can be decoded with ffdshow. Option -ffourcc MP4S will change it to MPS4 ("official" MS fourcc for mpeg4) and video will be playable on Windows by standard MS mpeg4 decoder. It may lead to side effects, if MS mpeg4 decoder is buggy. More info: http://www.itdp.de/mplayer-users/2005-03/msg00069.html

Other guides:
* http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-enc-images.html
* http://electron.mit.edu/~gsteele/ffmpeg/
* http://www.stillhq.com/extracted/howto-jpeg2mpeg/output.html

--------

==== Q: How to construct high-quality WebM movies? ====

A: {{wikipedia|WebM}} is the container for the {{wikipedia|VP8}} and {{wikipedia|VP9}} codecs. VP8 approaches H.264 in quality but differs from H.264 in that it is both patent + license free to use. Likewise, VP9 is an alternative to HEVC/H.265. The audio channel, if used, is encoded using the Ogg {{wikipedia|Vorbis}} codec for VP8 or Opus codec for VP9. Google has been promoting WebM as the standard fallback video format for HTML5 and is moving YouTube content to it. Firefox. Opera, and Google Chrome web browsers support WebM playback out of the box without need for extra plugins.

; Encoding:

First convert PNG or JPEG images to PPM format using the {{wikipedia|NetPBM}} tools.
: ''This step is pretty fast, I wouldn't worry about trying to parallelize it.''

<syntaxhighlight lang="bash">
for file in *.png ; do
pngtopnm $file > `basename $file .png`.ppm
done
</syntaxhighlight>
::or
<syntaxhighlight lang="bash">
for file in *.jpg ; do
jpegtopnm $file > `basename $file .jpg`.ppm
done
</syntaxhighlight>

Next create a Y4M stream containing a concatenation of the raw PPM frames:

<syntaxhighlight lang="bash">
cat *.ppm | ppmtoy4m -F 2:1 > all_frames.y4m
</syntaxhighlight>

* ''The "-F 2:1" option above sets the frame rate to 2 frames per second.''
* ''The <tt>ppmtoy4m</tt> program comes from the mjpegtools package.''
* ''For a very high number of frames you may wish to skip this step and create a named FIFO or stdin pipe directly into <tt>vpxenc</tt>, as the raw Y4M file can become extremely large.''

Finally, run the <tt>vpxenc</tt> encoder to create a high quality WebM video:
<syntaxhighlight lang="bash">
vpxenc all_frames.y4m -o my_movie.webm \
--passes=2 --threads=4 \
--best --target-bitrate=2000 --end-usage=vbr \
--auto-alt-ref=1 --verbose \
--minsection-pct=5 --maxsection-pct=800 \
--lag-in-frames=16 --kf-min-dist=0 --kf-max-dist=360 \
--static-thresh=0 --drop-frame=0 \
--min-q=0 --max-q=60
</syntaxhighlight>

* ''Set "--threads=" in the above command to the number of cores on your CPU. (or one less than the number of cores if you have many)''
* ''If you are encoding to VP8 (and not VP9) consider adding <tt>--token-parts=2</tt>''
* ''If <tt>vpxenc</tt> complains about "Invalid image format" either add "-S 420jpeg" to the <tt>ppmtoy4m</tt> options or add "--profile=1" to the <tt>vpxenc</tt> options.
* ''The <tt>vpxenc</tt> program comes from the vpx-tools package.''
* [http://www.webmproject.org/docs/encoder-parameters/#2-pass-best-quality-vbr-encoding Encoder parameters help page]

--------

==== Q: How to generate an animated GIF, FLI, or MNG? ====

A: for animations of 300 frames or less, animated GIF, [http://woodshole.er.usgs.gov/operations/modeling/flc.html FLI], or MNG formats are smaller files and better quality (i.e., frames stored as PNG not JPEG). Create with [http://www.lcdf.org/gifsicle/ gifsicle], [http://vento.pi.tu-berlin.de/fli.html ppm2fli], or for MNG use [http://www.imagemagick.org/script/convert.php ImageMagick's 'convert']. With more than 300 frames the players usually have memory issues. Animated GIFs are playable in any web browser of course and also in OpenOffice.

# with convert, many GIFs into one animated GIF:
convert -verbose -delay 20 -loop 0 snapshot*.png animation.gif

Notes: '-delay 20' means 20 hundreds of a second delay between each frame. The '-loop 0' flag lets it loop indefinitely. To loop only 3 times use '-loop 3' etc.

Sample script to export PNGs in GRASS (generates tmean_001.png .. tmean_708.png):

<syntaxhighlight lang=bash>
export i=0
for year in `seq 1950 2008` ; do
for m in `seq 1 12` ; do
i=`expr $i + 1`
# awk trick to generate wildcard-correct file names:
i=`echo $i | awk '{printf "%03d\n", $1}'`
r.out.png tmean_europe.$year.$m.avg out=tmean_$i.png
done
done
convert -verbose -delay 20 -loop 0 tmean_*.png tmean_animation.gif
</syntaxhighlight>

For presentations using a web browser, you can center the image on a blank white page, turn off any toolbars and go into full screen mode (F11 for Firefox or Chrome). Then Alt-Tab your way to the animation at the appropriate time and nobody knows you are using a web browser.

Example HTML for centering image:

<syntaxhighlight lang="html4strict">
<HTML>
<HEAD></HEAD>
<BODY>
 
<CENTER>
<IMG SRC="animation.gif">
</CENTER>
</BODY>
</HTML>
</syntaxhighlight>

--------

==== Q: How to create dynamic surface movies in NVIZ? ====

A: See slides from the FOSS4G 2006 workshop:
: http://skagit.meas.ncsu.edu/~helena/grasswork/foss4g/FOSS4G06WKSVisual4anim.odp

See also the example at the bottom of the [http://grass.osgeo.org/grass64/manuals/html64_user/nviz/nviz_panel_kanim.html NVIZ keyframe animator panel help page].

--------

==== Q: How to animate a time series of raster maps? ====

A: Use
* {{cmd|d.slide.show}} to view in a GRASS xmon (see {{cmd|d.mon}})
* {{cmd|r.out.mpeg}} to save to a MPEG-1 movie.
* GRASS GIS 6: {{cmd|xganim|version=64}} to view on screen in X-windows (with live start/stop/speed controls)
* GRASS GIS 7: {{cmd|g.gui.animation}} to view on screen in wxGUI (with live start/stop/speed controls)
: See also the wiki page for the [[wxGUI Animation Tool]]:
{{YouTube|Pi1ESk9iPwk|YouTube video introduction to the WxGUI Animation Tool}}

--------

==== Q: How to animate a series of maps, including decorations? ====

A: Use the xmon drivers ({{cmd|d.mon}})
to render the displays to a window via a shell script loop, then switch from the ''x0'' to ''PNG'' driver to save as a series of PNG or PPM images. Use a method answered in a FAQ above to combine those into a single animation.

--------

==== Q: How to create a screencast (for video tutorials etc)? ====

; GNU/Linux:
* Easist way: [http://recordmydesktop.sourceforge.net/ recordmydesktop] (user interfaces are 'qt-recordmydesktop' and 'gtk-recordmydesktop')
:* when uploading videos to Youtube created by <tt>recordmydesktop</tt> you need to encode the file eg. using <tt>mencoder</tt>
:- No sound, just video screencast, as AVI file:
::<pre>mencoder out.ogv -nosound -ovc lavc -o out.avi</pre>
:- With audio track, as AVI file:
::<pre>mencoder out.ogv -ovc xvid -oac mp3lame -xvidencopts pass=1 -o out.avi</pre>
:- With audio track, as WebM file:
::<pre>ffmpeg -i out.ogv -b:a 24k out.webm</pre>
:- With audio track, as MP4 file:
::<pre>ffmpeg -i out.ogv -ar 44100 -ab 96k -coder ac -strict experimental -vcodec libx264 out.mp4</pre>
* [https://wiki.ubuntu.com/MeetingLogs/openweekfeisty/screencast Ubuntu discussion over screencast]
* [https://wiki.ubuntu.com/ScreencastTeam/RecordingScreencasts Ubuntu screencast how to]
* [http://www.unixuser.org/~euske/vnc2swf/ vnc2swf] to record a flash movie ([http://www-pool.math.tu-berlin.de/~soeren/grass/modules/screenshots/vnc2swf_usage.html short vnc2swf usage tutorial])
* [http://www.maartenbaert.be/simplescreenrecorder/ SimpleScreenRecorder]

; MacOSX:
* [http://crschmidt.net/blog/archives/203/openlayers-screencast/ Chris Schmidt experience on MacOSX]

; MS Windows:
* [http://camstudio.org/ Camstudio]

-----

==== Examples ====
''Using the [http://grass.osgeo.org/download/data6.php Spearfish] sample dataset''

===== Animate on screen =====
''Using a shell script''

<syntaxhighlight lang=bash>
GRASS_WIDTH=640 GRASS_HEIGHT=500 \
d.mon start=x0

g.region rast=elevation.dem n=4927830 s=4912980

for NUM in `seq 1 0.5 10` ; do
d.rast elevation.dem --quiet
d.vect roads
d.vect bugsites icon=basic/circle col=black fcol=green size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph

echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
done
</syntaxhighlight>

===== Animate to files =====
''Using [http://grass.osgeo.org/grass64/manuals/html64_user/pngdriver.html GRASS's PNG driver]''

<syntaxhighlight lang=bash>
# set up base map
g.region rast=elevation.dem n=4927830 s=4912980
export GRASS_WIDTH=640
export GRASS_HEIGHT=500
export GRASS_PNGFILE=base.ppm
export GRASS_PNG_READ=FALSE
d.mon start=PNG
d.rast elevation.dem
d.vect roads
echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
d.mon stop=PNG

# loop to draw variable parts as new frames
i=0
export GRASS_PNG_READ=TRUE
for NUM in `seq 1 0.5 10` ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
echo "Processing frame $FRAMENUM ..."
export GRASS_PNGFILE="ganim_${FRAMENUM}.ppm"
cp base.ppm "$GRASS_PNGFILE"
d.mon start=PNG --quiet
d.vect bugsites icon=basic/circle col=black fcol=green \
size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph
d.mon stop=PNG --quiet
done

# repeat the last frame a few times
for EXTRA in 1 2 3 ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
cp base.ppm "ganim_${FRAMENUM}.ppm"
done
\rm base.ppm
</syntaxhighlight>

===== Encode as animated GIF =====

* ''Using the [http://www.lcdf.org/gifsicle/ gifsicle] encoder''
* ''View in a web browser or with ImageMagick's [http://www.imagemagick.org/Usage/basics/#animate animate]''
* see [[Movies#Q:_How_to_generate_an_animated_GIF.2C_FLI.2C_or_MNG.3F]] above

'''A) Using gifsicle:'''

Advantage: the delay is coded in the animated GIF file and rendered properly e.g. in a Web browser.

<syntaxhighlight lang=bash>
# Imagemagick's `convert` program may do a better job at
# 256 color quantization than NetPBM's <tt>ppmquant</tt>.
for IMG in ganim*.ppm ; do # convert each frame from PPM to GIF
#ppmquant 256 $IMG | ppmtogif > "`basename $IMG .ppm`.gif"
convert "$IMG" "`basename $IMG .ppm`.gif"
done

gifsicle -O2 --delay 20 --no-loopcount --colors 256 ganim_*.gif > ganim.gif

\rm ganim_0*.gif
</syntaxhighlight>

'''B) Using ImageMagick's convert and [http://www.imagemagick.org/Usage/basics/#animate animate]:'''
<syntaxhighlight lang=bash>
# Generate "film" (create animated GIF, 750 millisec between frames):
convert -delay 75 lst_monthly_avg_film.00*.png lst_monthly_avg_2002_2012_zoom_film.gif

# play animation (delay 500 millisec between frames, pause for 5 seconds at the end)
animate -delay 50 -pause 5 lst_monthly_avg_2002_2012_zoom_film.gif
</syntaxhighlight>

===== Encode as animated PNG (MNG) =====
* ''Using the [http://www.imagemagick.org ImageMagick] "convert" encoder''
* ''View with ImageMagick's "display"''

<pre>
convert -delay 20 ganim_*.ppm ganim.mng
</pre>

===== Encode as FLI =====
* ''Using the [http://vento.pi.tu-berlin.de/fli.html ppm2fli] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or [http://xanim.polter.net xanim]''

<syntaxhighlight lang=bash>
ls ganim_*.ppm > frames.txt
ppm2fli -g"${GRASS_WIDTH}x${GRASS_HEIGHT}" -O -s 15 frames.txt ganim.fli
\rm frames.txt
</syntaxhighlight>

===== Encode as Flash =====

* ''Using the [http://www.swftools.org/ SWFTools] <tt>png2swf</tt> encoder''
* ''View with [http://www.mplayerhq.hu mplayer]''

png2swf -o outfile.swf *.png -r 1

===== Encode as MPEG-4 Xvid =====
* ''Using the [http://www.mplayerhq.hu mencoder] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or similar''

<syntaxhighlight lang=bash>
# doesn't like GRASS's PPMs, convert to PNG
for IMG in ganim*.ppm ; do
pnmtopng $IMG > "`basename $IMG .ppm`.png"
done
mencoder "mf://ganim_*.png" -mf "type=png:fps=5" -o ganim.avi \
-ovc "xvid" -xvidencopts "bitrate=1024"
\rm ganim_0*.png
</syntaxhighlight>

==== Publish your cool videos in this Wiki ====
We have installed a Mediawiki widget extension and created a special template for including Youtube videos ([[Template:YouTube]]). See also [[WxGUI_Modeler#Video_tutorials|wxGUI Modeler]] page for the usage.

The overall goal is to have many as possible '''[http://www.youtube.com/results?search_query=grass+gis GRASS GIS videos on YouTube]'''

[[Category:FAQ]]
[[Category:Visualization]]

Movies

2025-11-04T03:06:33Z

Hamish: /* Encode as animated GIF */

==== Q: How to construct high-quality MPEG-4 movies (animations) from a series of still frames? ====

''Maris wrote'':

Encode all .png files in directory to out.avi with 15 frames per second:
<pre>
mencoder "mf:///path/to/files/*.png" -mf fps=15 -o out.avi # put -ovc here, see next lines
</pre>
<pre>
## for DivX - libavcodec MPEG 4 (DivX5), use:
-ovc lavc -lavcopts vcodec=mpeg4:vhq:vbitrate=1800 -ffourcc MP4S
## for XviD, use:
-ovc xvid -xvidencopts bitrate=1024 -xvidencopts pass=2 -xvidencopts cartoon
## for DivX4, use:
-ovc divx4 -divx4opts q=5
## for x264 use:
-ovc x264 -x264encopts pass=3 -x264encopts preset=veryslow -x264encopts tune=animation
</pre>

Uses mplayer's encoder. Choose one of encoding formats and append to end
of "mencoder" line. Unfortunately only way to know quality settings is by
encoding, watching result and re-encoding with different params.

Notice - by default encoding with lavc will set video fourcc to FMP4, that can be decoded with ffdshow. Option -ffourcc MP4S will change it to MPS4 ("official" MS fourcc for mpeg4) and video will be playable on Windows by standard MS mpeg4 decoder. It may lead to side effects, if MS mpeg4 decoder is buggy. More info: http://www.itdp.de/mplayer-users/2005-03/msg00069.html

Other guides:
* http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-enc-images.html
* http://electron.mit.edu/~gsteele/ffmpeg/
* http://www.stillhq.com/extracted/howto-jpeg2mpeg/output.html

--------

==== Q: How to construct high-quality WebM movies? ====

A: {{wikipedia|WebM}} is the container for the {{wikipedia|VP8}} and {{wikipedia|VP9}} codecs. VP8 approaches H.264 in quality but differs from H.264 in that it is both patent + license free to use. Likewise, VP9 is an alternative to HEVC/H.265. The audio channel, if used, is encoded using the Ogg {{wikipedia|Vorbis}} codec for VP8 or Opus codec for VP9. Google has been promoting WebM as the standard fallback video format for HTML5 and is moving YouTube content to it. Firefox. Opera, and Google Chrome web browsers support WebM playback out of the box without need for extra plugins.

; Encoding:

First convert PNG or JPEG images to PPM format using the {{wikipedia|NetPBM}} tools.
: ''This step is pretty fast, I wouldn't worry about trying to parallelize it.''

<source lang="bash">
for file in *.png ; do
pngtopnm $file > `basename $file .png`.ppm
done
</source>
::or
<source lang="bash">
for file in *.jpg ; do
jpegtopnm $file > `basename $file .jpg`.ppm
done
</source>

Next create a Y4M stream containing a concatenation of the raw PPM frames:

<source lang="bash">
cat *.ppm | ppmtoy4m -F 2:1 > all_frames.y4m
</source>

* ''The "-F 2:1" option above sets the frame rate to 2 frames per second.''
* ''The <tt>ppmtoy4m</tt> program comes from the mjpegtools package.''
* ''For a very high number of frames you may wish to skip this step and create a named FIFO or stdin pipe directly into <tt>vpxenc</tt>, as the raw Y4M file can become extremely large.''

Finally, run the <tt>vpxenc</tt> encoder to create a high quality WebM video:
<source lang="bash">
vpxenc all_frames.y4m -o my_movie.webm \
--passes=2 --threads=4 \
--best --target-bitrate=2000 --end-usage=vbr \
--auto-alt-ref=1 --verbose \
--minsection-pct=5 --maxsection-pct=800 \
--lag-in-frames=16 --kf-min-dist=0 --kf-max-dist=360 \
--static-thresh=0 --drop-frame=0 \
--min-q=0 --max-q=60
</source>

* ''Set "--threads=" in the above command to the number of cores on your CPU. (or one less than the number of cores if you have many)''
* ''If you are encoding to VP8 (and not VP9) consider adding <tt>--token-parts=2</tt>''
* ''If <tt>vpxenc</tt> complains about "Invalid image format" either add "-S 420jpeg" to the <tt>ppmtoy4m</tt> options or add "--profile=1" to the <tt>vpxenc</tt> options.
* ''The <tt>vpxenc</tt> program comes from the vpx-tools package.''
* [http://www.webmproject.org/docs/encoder-parameters/#2-pass-best-quality-vbr-encoding Encoder parameters help page]

--------

==== Q: How to generate an animated GIF, FLI, or MNG? ====

A: for animations of 300 frames or less, animated GIF, [http://woodshole.er.usgs.gov/operations/modeling/flc.html FLI], or MNG formats are smaller files and better quality (i.e., frames stored as PNG not JPEG). Create with [http://www.lcdf.org/gifsicle/ gifsicle], [http://vento.pi.tu-berlin.de/fli.html ppm2fli], or for MNG use [http://www.imagemagick.org/script/convert.php ImageMagick's 'convert']. With more than 300 frames the players usually have memory issues. Animated GIFs are playable in any web browser of course and also in OpenOffice.

# with convert, many GIFs into one animated GIF:
convert -verbose -delay 20 -loop 0 snapshot*.png animation.gif

Notes: '-delay 20' means 20 hundreds of a second delay between each frame. The '-loop 0' flag lets it loop indefinitely. To loop only 3 times use '-loop 3' etc.

Sample script to export PNGs in GRASS (generates tmean_001.png .. tmean_708.png):

<source lang=bash>
export i=0
for year in `seq 1950 2008` ; do
for m in `seq 1 12` ; do
i=`expr $i + 1`
# awk trick to generate wildcard-correct file names:
i=`echo $i | awk '{printf "%03d\n", $1}'`
r.out.png tmean_europe.$year.$m.avg out=tmean_$i.png
done
done
convert -verbose -delay 20 -loop 0 tmean_*.png tmean_animation.gif
</source>

For presentations using a web browser, you can center the image on a blank white page, turn off any toolbars and go into full screen mode (F11 for Firefox or Chrome). Then Alt-Tab your way to the animation at the appropriate time and nobody knows you are using a web browser.

Example HTML for centering image:

<source lang="html4strict">
<HTML>
<HEAD></HEAD>
<BODY>
 
<CENTER>
<IMG SRC="animation.gif">
</CENTER>
</BODY>
</HTML>
</source>

--------

==== Q: How to create dynamic surface movies in NVIZ? ====

A: See slides from the FOSS4G 2006 workshop:
: http://skagit.meas.ncsu.edu/~helena/grasswork/foss4g/FOSS4G06WKSVisual4anim.odp

See also the example at the bottom of the [http://grass.osgeo.org/grass64/manuals/html64_user/nviz/nviz_panel_kanim.html NVIZ keyframe animator panel help page].

--------

==== Q: How to animate a time series of raster maps? ====

A: Use
* {{cmd|d.slide.show}} to view in a GRASS xmon (see {{cmd|d.mon}})
* {{cmd|r.out.mpeg}} to save to a MPEG-1 movie.
* GRASS GIS 6: {{cmd|xganim|version=64}} to view on screen in X-windows (with live start/stop/speed controls)
* GRASS GIS 7: {{cmd|g.gui.animation}} to view on screen in wxGUI (with live start/stop/speed controls)
: See also the wiki page for the [[wxGUI Animation Tool]]:
{{YouTube|Pi1ESk9iPwk|YouTube video introduction to the WxGUI Animation Tool}}

--------

==== Q: How to animate a series of maps, including decorations? ====

A: Use the xmon drivers ({{cmd|d.mon}})
to render the displays to a window via a shell script loop, then switch from the ''x0'' to ''PNG'' driver to save as a series of PNG or PPM images. Use a method answered in a FAQ above to combine those into a single animation.

--------

==== Q: How to create a screencast (for video tutorials etc)? ====

; GNU/Linux:
* Easist way: [http://recordmydesktop.sourceforge.net/ recordmydesktop] (user interfaces are 'qt-recordmydesktop' and 'gtk-recordmydesktop')
:* when uploading videos to Youtube created by <tt>recordmydesktop</tt> you need to encode the file eg. using <tt>mencoder</tt>
:- No sound, just video screencast, as AVI file:
::<pre>mencoder out.ogv -nosound -ovc lavc -o out.avi</pre>
:- With audio track, as AVI file:
::<pre>mencoder out.ogv -ovc xvid -oac mp3lame -xvidencopts pass=1 -o out.avi</pre>
:- With audio track, as WebM file:
::<pre>ffmpeg -i out.ogv -b:a 24k out.webm</pre>
:- With audio track, as MP4 file:
::<pre>ffmpeg -i out.ogv -ar 44100 -ab 96k -coder ac -strict experimental -vcodec libx264 out.mp4</pre>
* [https://wiki.ubuntu.com/MeetingLogs/openweekfeisty/screencast Ubuntu discussion over screencast]
* [https://wiki.ubuntu.com/ScreencastTeam/RecordingScreencasts Ubuntu screencast how to]
* [http://www.unixuser.org/~euske/vnc2swf/ vnc2swf] to record a flash movie ([http://www-pool.math.tu-berlin.de/~soeren/grass/modules/screenshots/vnc2swf_usage.html short vnc2swf usage tutorial])
* [http://www.maartenbaert.be/simplescreenrecorder/ SimpleScreenRecorder]

; MacOSX:
* [http://crschmidt.net/blog/archives/203/openlayers-screencast/ Chris Schmidt experience on MacOSX]

; MS Windows:
* [http://camstudio.org/ Camstudio]

-----

==== Examples ====
''Using the [http://grass.osgeo.org/download/data6.php Spearfish] sample dataset''

===== Animate on screen =====
''Using a shell script''

<source lang=bash>
GRASS_WIDTH=640 GRASS_HEIGHT=500 \
d.mon start=x0

g.region rast=elevation.dem n=4927830 s=4912980

for NUM in `seq 1 0.5 10` ; do
d.rast elevation.dem --quiet
d.vect roads
d.vect bugsites icon=basic/circle col=black fcol=green size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph

echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
done
</source>

===== Animate to files =====
''Using [http://grass.osgeo.org/grass64/manuals/html64_user/pngdriver.html GRASS's PNG driver]''

<source lang=bash>
# set up base map
g.region rast=elevation.dem n=4927830 s=4912980
export GRASS_WIDTH=640
export GRASS_HEIGHT=500
export GRASS_PNGFILE=base.ppm
export GRASS_PNG_READ=FALSE
d.mon start=PNG
d.rast elevation.dem
d.vect roads
echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
d.mon stop=PNG

# loop to draw variable parts as new frames
i=0
export GRASS_PNG_READ=TRUE
for NUM in `seq 1 0.5 10` ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
echo "Processing frame $FRAMENUM ..."
export GRASS_PNGFILE="ganim_${FRAMENUM}.ppm"
cp base.ppm "$GRASS_PNGFILE"
d.mon start=PNG --quiet
d.vect bugsites icon=basic/circle col=black fcol=green \
size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph
d.mon stop=PNG --quiet
done

# repeat the last frame a few times
for EXTRA in 1 2 3 ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
cp base.ppm "ganim_${FRAMENUM}.ppm"
done
\rm base.ppm
</source>

===== Encode as animated GIF =====

* ''Using the [http://www.lcdf.org/gifsicle/ gifsicle] encoder''
* ''View in a web browser or with ImageMagick's [http://www.imagemagick.org/Usage/basics/#animate animate]''
* see [[Movies#Q:_How_to_generate_an_animated_GIF.2C_FLI.2C_or_MNG.3F]] above

'''A) Using gifsicle:'''

Advantage: the delay is coded in the animated GIF file and rendered properly e.g. in a Web browser.

<source lang=bash>
# Imagemagick's `convert` program may do a better job at
# 256 color quantization than NetPBM's <tt>ppmquant</tt>.
for IMG in ganim*.ppm ; do # convert each frame from PPM to GIF
#ppmquant 256 $IMG | ppmtogif > "`basename $IMG .ppm`.gif"
convert "$IMG" "`basename $IMG .ppm`.gif"
done

gifsicle -O2 --delay 20 --no-loopcount --colors 256 ganim_*.gif > ganim.gif

\rm ganim_0*.gif
</source>

'''B) Using ImageMagick's convert and [http://www.imagemagick.org/Usage/basics/#animate animate]:'''
<source lang=bash>
# Generate "film" (create animated GIF, 750 millisec between frames):
convert -delay 75 lst_monthly_avg_film.00*.png lst_monthly_avg_2002_2012_zoom_film.gif

# play animation (delay 500 millisec between frames, pause for 5 seconds at the end)
animate -delay 50 -pause 5 lst_monthly_avg_2002_2012_zoom_film.gif
</source>

===== Encode as animated PNG (MNG) =====
* ''Using the [http://www.imagemagick.org ImageMagick] "convert" encoder''
* ''View with ImageMagick's "display"''

<pre>
convert -delay 20 ganim_*.ppm ganim.mng
</pre>

===== Encode as FLI =====
* ''Using the [http://vento.pi.tu-berlin.de/fli.html ppm2fli] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or [http://xanim.polter.net xanim]''

<source lang=bash>
ls ganim_*.ppm > frames.txt
ppm2fli -g"${GRASS_WIDTH}x${GRASS_HEIGHT}" -O -s 15 frames.txt ganim.fli
\rm frames.txt
</source>

===== Encode as Flash =====

* ''Using the [http://www.swftools.org/ SWFTools] <tt>png2swf</tt> encoder''
* ''View with [http://www.mplayerhq.hu mplayer]''

png2swf -o outfile.swf *.png -r 1

===== Encode as MPEG-4 Xvid =====
* ''Using the [http://www.mplayerhq.hu mencoder] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or similar''

<source lang=bash>
# doesn't like GRASS's PPMs, convert to PNG
for IMG in ganim*.ppm ; do
pnmtopng $IMG > "`basename $IMG .ppm`.png"
done
mencoder "mf://ganim_*.png" -mf "type=png:fps=5" -o ganim.avi \
-ovc "xvid" -xvidencopts "bitrate=1024"
\rm ganim_0*.png
</source>

==== Publish your cool videos in this Wiki ====
We have installed a Mediawiki widget extension and created a special template for including Youtube videos ([[Template:YouTube]]). See also [[WxGUI_Modeler#Video_tutorials|wxGUI Modeler]] page for the usage.

The overall goal is to have many as possible '''[http://www.youtube.com/results?search_query=grass+gis GRASS GIS videos on YouTube]'''

[[Category:FAQ]]
[[Category:Visualization]]

Movies

2025-11-04T03:04:42Z

Hamish: /* Encode as animated GIF */

==== Q: How to construct high-quality MPEG-4 movies (animations) from a series of still frames? ====

''Maris wrote'':

Encode all .png files in directory to out.avi with 15 frames per second:
<pre>
mencoder "mf:///path/to/files/*.png" -mf fps=15 -o out.avi # put -ovc here, see next lines
</pre>
<pre>
## for DivX - libavcodec MPEG 4 (DivX5), use:
-ovc lavc -lavcopts vcodec=mpeg4:vhq:vbitrate=1800 -ffourcc MP4S
## for XviD, use:
-ovc xvid -xvidencopts bitrate=1024 -xvidencopts pass=2 -xvidencopts cartoon
## for DivX4, use:
-ovc divx4 -divx4opts q=5
## for x264 use:
-ovc x264 -x264encopts pass=3 -x264encopts preset=veryslow -x264encopts tune=animation
</pre>

Uses mplayer's encoder. Choose one of encoding formats and append to end
of "mencoder" line. Unfortunately only way to know quality settings is by
encoding, watching result and re-encoding with different params.

Notice - by default encoding with lavc will set video fourcc to FMP4, that can be decoded with ffdshow. Option -ffourcc MP4S will change it to MPS4 ("official" MS fourcc for mpeg4) and video will be playable on Windows by standard MS mpeg4 decoder. It may lead to side effects, if MS mpeg4 decoder is buggy. More info: http://www.itdp.de/mplayer-users/2005-03/msg00069.html

Other guides:
* http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-enc-images.html
* http://electron.mit.edu/~gsteele/ffmpeg/
* http://www.stillhq.com/extracted/howto-jpeg2mpeg/output.html

--------

==== Q: How to construct high-quality WebM movies? ====

A: {{wikipedia|WebM}} is the container for the {{wikipedia|VP8}} and {{wikipedia|VP9}} codecs. VP8 approaches H.264 in quality but differs from H.264 in that it is both patent + license free to use. Likewise, VP9 is an alternative to HEVC/H.265. The audio channel, if used, is encoded using the Ogg {{wikipedia|Vorbis}} codec for VP8 or Opus codec for VP9. Google has been promoting WebM as the standard fallback video format for HTML5 and is moving YouTube content to it. Firefox. Opera, and Google Chrome web browsers support WebM playback out of the box without need for extra plugins.

; Encoding:

First convert PNG or JPEG images to PPM format using the {{wikipedia|NetPBM}} tools.
: ''This step is pretty fast, I wouldn't worry about trying to parallelize it.''

<source lang="bash">
for file in *.png ; do
pngtopnm $file > `basename $file .png`.ppm
done
</source>
::or
<source lang="bash">
for file in *.jpg ; do
jpegtopnm $file > `basename $file .jpg`.ppm
done
</source>

Next create a Y4M stream containing a concatenation of the raw PPM frames:

<source lang="bash">
cat *.ppm | ppmtoy4m -F 2:1 > all_frames.y4m
</source>

* ''The "-F 2:1" option above sets the frame rate to 2 frames per second.''
* ''The <tt>ppmtoy4m</tt> program comes from the mjpegtools package.''
* ''For a very high number of frames you may wish to skip this step and create a named FIFO or stdin pipe directly into <tt>vpxenc</tt>, as the raw Y4M file can become extremely large.''

Finally, run the <tt>vpxenc</tt> encoder to create a high quality WebM video:
<source lang="bash">
vpxenc all_frames.y4m -o my_movie.webm \
--passes=2 --threads=4 \
--best --target-bitrate=2000 --end-usage=vbr \
--auto-alt-ref=1 --verbose \
--minsection-pct=5 --maxsection-pct=800 \
--lag-in-frames=16 --kf-min-dist=0 --kf-max-dist=360 \
--static-thresh=0 --drop-frame=0 \
--min-q=0 --max-q=60
</source>

* ''Set "--threads=" in the above command to the number of cores on your CPU. (or one less than the number of cores if you have many)''
* ''If you are encoding to VP8 (and not VP9) consider adding <tt>--token-parts=2</tt>''
* ''If <tt>vpxenc</tt> complains about "Invalid image format" either add "-S 420jpeg" to the <tt>ppmtoy4m</tt> options or add "--profile=1" to the <tt>vpxenc</tt> options.
* ''The <tt>vpxenc</tt> program comes from the vpx-tools package.''
* [http://www.webmproject.org/docs/encoder-parameters/#2-pass-best-quality-vbr-encoding Encoder parameters help page]

--------

==== Q: How to generate an animated GIF, FLI, or MNG? ====

A: for animations of 300 frames or less, animated GIF, [http://woodshole.er.usgs.gov/operations/modeling/flc.html FLI], or MNG formats are smaller files and better quality (i.e., frames stored as PNG not JPEG). Create with [http://www.lcdf.org/gifsicle/ gifsicle], [http://vento.pi.tu-berlin.de/fli.html ppm2fli], or for MNG use [http://www.imagemagick.org/script/convert.php ImageMagick's 'convert']. With more than 300 frames the players usually have memory issues. Animated GIFs are playable in any web browser of course and also in OpenOffice.

# with convert, many GIFs into one animated GIF:
convert -verbose -delay 20 -loop 0 snapshot*.png animation.gif

Notes: '-delay 20' means 20 hundreds of a second delay between each frame. The '-loop 0' flag lets it loop indefinitely. To loop only 3 times use '-loop 3' etc.

Sample script to export PNGs in GRASS (generates tmean_001.png .. tmean_708.png):

<source lang=bash>
export i=0
for year in `seq 1950 2008` ; do
for m in `seq 1 12` ; do
i=`expr $i + 1`
# awk trick to generate wildcard-correct file names:
i=`echo $i | awk '{printf "%03d\n", $1}'`
r.out.png tmean_europe.$year.$m.avg out=tmean_$i.png
done
done
convert -verbose -delay 20 -loop 0 tmean_*.png tmean_animation.gif
</source>

For presentations using a web browser, you can center the image on a blank white page, turn off any toolbars and go into full screen mode (F11 for Firefox or Chrome). Then Alt-Tab your way to the animation at the appropriate time and nobody knows you are using a web browser.

Example HTML for centering image:

<source lang="html4strict">
<HTML>
<HEAD></HEAD>
<BODY>
 
<CENTER>
<IMG SRC="animation.gif">
</CENTER>
</BODY>
</HTML>
</source>

--------

==== Q: How to create dynamic surface movies in NVIZ? ====

A: See slides from the FOSS4G 2006 workshop:
: http://skagit.meas.ncsu.edu/~helena/grasswork/foss4g/FOSS4G06WKSVisual4anim.odp

See also the example at the bottom of the [http://grass.osgeo.org/grass64/manuals/html64_user/nviz/nviz_panel_kanim.html NVIZ keyframe animator panel help page].

--------

==== Q: How to animate a time series of raster maps? ====

A: Use
* {{cmd|d.slide.show}} to view in a GRASS xmon (see {{cmd|d.mon}})
* {{cmd|r.out.mpeg}} to save to a MPEG-1 movie.
* GRASS GIS 6: {{cmd|xganim|version=64}} to view on screen in X-windows (with live start/stop/speed controls)
* GRASS GIS 7: {{cmd|g.gui.animation}} to view on screen in wxGUI (with live start/stop/speed controls)
: See also the wiki page for the [[wxGUI Animation Tool]]:
{{YouTube|Pi1ESk9iPwk|YouTube video introduction to the WxGUI Animation Tool}}

--------

==== Q: How to animate a series of maps, including decorations? ====

A: Use the xmon drivers ({{cmd|d.mon}})
to render the displays to a window via a shell script loop, then switch from the ''x0'' to ''PNG'' driver to save as a series of PNG or PPM images. Use a method answered in a FAQ above to combine those into a single animation.

--------

==== Q: How to create a screencast (for video tutorials etc)? ====

; GNU/Linux:
* Easist way: [http://recordmydesktop.sourceforge.net/ recordmydesktop] (user interfaces are 'qt-recordmydesktop' and 'gtk-recordmydesktop')
:* when uploading videos to Youtube created by <tt>recordmydesktop</tt> you need to encode the file eg. using <tt>mencoder</tt>
:- No sound, just video screencast, as AVI file:
::<pre>mencoder out.ogv -nosound -ovc lavc -o out.avi</pre>
:- With audio track, as AVI file:
::<pre>mencoder out.ogv -ovc xvid -oac mp3lame -xvidencopts pass=1 -o out.avi</pre>
:- With audio track, as WebM file:
::<pre>ffmpeg -i out.ogv -b:a 24k out.webm</pre>
:- With audio track, as MP4 file:
::<pre>ffmpeg -i out.ogv -ar 44100 -ab 96k -coder ac -strict experimental -vcodec libx264 out.mp4</pre>
* [https://wiki.ubuntu.com/MeetingLogs/openweekfeisty/screencast Ubuntu discussion over screencast]
* [https://wiki.ubuntu.com/ScreencastTeam/RecordingScreencasts Ubuntu screencast how to]
* [http://www.unixuser.org/~euske/vnc2swf/ vnc2swf] to record a flash movie ([http://www-pool.math.tu-berlin.de/~soeren/grass/modules/screenshots/vnc2swf_usage.html short vnc2swf usage tutorial])
* [http://www.maartenbaert.be/simplescreenrecorder/ SimpleScreenRecorder]

; MacOSX:
* [http://crschmidt.net/blog/archives/203/openlayers-screencast/ Chris Schmidt experience on MacOSX]

; MS Windows:
* [http://camstudio.org/ Camstudio]

-----

==== Examples ====
''Using the [http://grass.osgeo.org/download/data6.php Spearfish] sample dataset''

===== Animate on screen =====
''Using a shell script''

<source lang=bash>
GRASS_WIDTH=640 GRASS_HEIGHT=500 \
d.mon start=x0

g.region rast=elevation.dem n=4927830 s=4912980

for NUM in `seq 1 0.5 10` ; do
d.rast elevation.dem --quiet
d.vect roads
d.vect bugsites icon=basic/circle col=black fcol=green size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph

echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
done
</source>

===== Animate to files =====
''Using [http://grass.osgeo.org/grass64/manuals/html64_user/pngdriver.html GRASS's PNG driver]''

<source lang=bash>
# set up base map
g.region rast=elevation.dem n=4927830 s=4912980
export GRASS_WIDTH=640
export GRASS_HEIGHT=500
export GRASS_PNGFILE=base.ppm
export GRASS_PNG_READ=FALSE
d.mon start=PNG
d.rast elevation.dem
d.vect roads
echo "Spearfish, South Dakota" | d.text size=4.5 at=1,2.25
d.barscale at=59.8,93.6 tcol=grey
d.mon stop=PNG

# loop to draw variable parts as new frames
i=0
export GRASS_PNG_READ=TRUE
for NUM in `seq 1 0.5 10` ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
echo "Processing frame $FRAMENUM ..."
export GRASS_PNGFILE="ganim_${FRAMENUM}.ppm"
cp base.ppm "$GRASS_PNGFILE"
d.mon start=PNG --quiet
d.vect bugsites icon=basic/circle col=black fcol=green \
size=`echo "21-($NUM * 2)" | bc`

echo "symbol basic/arrow2 25 80 `echo "10*$NUM" | bc` 80:80:80 125:125:255" | d.graph
d.mon stop=PNG --quiet
done

# repeat the last frame a few times
for EXTRA in 1 2 3 ; do
i=`expr $i + 1`
FRAMENUM=`echo $i | awk '{printf("%03d", $1)}'`
cp base.ppm "ganim_${FRAMENUM}.ppm"
done
\rm base.ppm
</source>

===== Encode as animated GIF =====

* ''Using the [http://www.lcdf.org/gifsicle/ gifsicle] encoder''
* ''View in a web browser or with ImageMagick's [http://www.imagemagick.org/Usage/basics/#animate animate]''
* see [[Movies#Q:_How_to_generate_an_animated_GIF.2C_FLI.2C_or_MNG.3F]] above

'''A) Using gifsicle:'''

Advantage: the delay is coded in the animated GIF file and rendered properly e.g. in a Web browser.

<source lang=bash>
# Imagemagick's `convert` program may do a better job at
# 256 color quantization than NetPBM's <tt>ppmquant</tt>.
for IMG in ganim*.ppm ; do # convert each frame from PPM to GIF
#ppmquant 256 $IMG | ppmtogif > "`basename $IMG .ppm`.gif"
convert "$IMG" "`basename $IMG .ppm`.gif"
done

gifsicle -O2 --delay 20 --no-loopcount --colors 256 ganim_*.gif > ganim.gif

\rm ganim_0*.gif
</source>

'''B) Using ImageMagick's convert and [http://www.imagemagick.org/Usage/basics/#animate animate]:'''
<source lang=bash>
# Generate "film" (create animated GIF):
convert lst_monthly_avg_film.00*.png lst_monthly_avg_2002_2012_zoom_film.gif

# play animation (delay 50millisec between frames, pause for 5 seconds at the end)
animate -delay 50 -pause 5 lst_monthly_avg_2002_2012_zoom_film.gif
</source>

===== Encode as animated PNG (MNG) =====
* ''Using the [http://www.imagemagick.org ImageMagick] "convert" encoder''
* ''View with ImageMagick's "display"''

<pre>
convert -delay 20 ganim_*.ppm ganim.mng
</pre>

===== Encode as FLI =====
* ''Using the [http://vento.pi.tu-berlin.de/fli.html ppm2fli] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or [http://xanim.polter.net xanim]''

<source lang=bash>
ls ganim_*.ppm > frames.txt
ppm2fli -g"${GRASS_WIDTH}x${GRASS_HEIGHT}" -O -s 15 frames.txt ganim.fli
\rm frames.txt
</source>

===== Encode as Flash =====

* ''Using the [http://www.swftools.org/ SWFTools] <tt>png2swf</tt> encoder''
* ''View with [http://www.mplayerhq.hu mplayer]''

png2swf -o outfile.swf *.png -r 1

===== Encode as MPEG-4 Xvid =====
* ''Using the [http://www.mplayerhq.hu mencoder] encoder''
* ''View with [http://www.mplayerhq.hu mplayer] or similar''

<source lang=bash>
# doesn't like GRASS's PPMs, convert to PNG
for IMG in ganim*.ppm ; do
pnmtopng $IMG > "`basename $IMG .ppm`.png"
done
mencoder "mf://ganim_*.png" -mf "type=png:fps=5" -o ganim.avi \
-ovc "xvid" -xvidencopts "bitrate=1024"
\rm ganim_0*.png
</source>

==== Publish your cool videos in this Wiki ====
We have installed a Mediawiki widget extension and created a special template for including Youtube videos ([[Template:YouTube]]). See also [[WxGUI_Modeler#Video_tutorials|wxGUI Modeler]] page for the usage.

The overall goal is to have many as possible '''[http://www.youtube.com/results?search_query=grass+gis GRASS GIS videos on YouTube]'''

[[Category:FAQ]]
[[Category:Visualization]]

R.sun

2025-10-29T05:12:52Z

Hamish: /* Linke turbidity factor values */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. Note that this may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Linke turbidity factor values ==

* See Remund et al. (2003) appendix/table 1 for a list of monthly values for many sites around the world:
: https://www.researchgate.net/publication/42358500_Worldwide_Linke_turbidity_information

* A good description and recipe for calculating it here: (section 3.3.4, page 22)
: https://atmosphere.copernicus.eu/sites/default/files/2022-01/CAMS2_73_2021SC1_D3.2.1_2021_UserGuide_v1.pdf

* Download 5' global dataset for all months of 2010 as LL/WGS84 GeoTIFF here:
: https://www.soda-pro.com/help/general-knowledge/linke-turbidity-factor
:: '''''"The Linke turbidity factor (TL, for an air mass equal to 2) is a very convenient approximation to model the atmospheric absorption and scattering of the solar radiation under clear skies. It describes the optical thickness of the atmosphere due to both the absorption by the water vapor and the absorption and scattering by the aerosol particles relative to a dry and clean atmosphere. It summarizes the turbidity of the atmosphere, and hence the attenuation of the direct beam solar radiation. The larger TL, the larger the attenuation of the radiation by the clear sky atmosphere.'''''

:: '''''The TL denotes the transparency of the cloudless atmosphere. If the sky were perfectly dry and clean, TL would be equal to 1. When the sky is deep blue, the TL is just above 1 and still very small. In summer, in Europe, the water vapour is often large and the blue sky is close to white -- the TL is larger than 3. In turbid atmosphere, e.g. in polluted cities, the TL is close to 6 - 7."'''''

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get an estimated curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-10-29T05:11:09Z

Hamish: /* Linke turbidity factor values */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. Note that this may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Linke turbidity factor values ==

* See Remund et al. (2003) appendix/table 1 for a list of monthly values for many sites around the world:
: https://www.researchgate.net/publication/42358500_Worldwide_Linke_turbidity_information

* A good description and recipe for calculating it here: (section 3.3.4, page 22)
: https://atmosphere.copernicus.eu/sites/default/files/2022-01/CAMS2_73_2021SC1_D3.2.1_2021_UserGuide_v1.pdf

* Download 5' global dataset for all months of 2010 as LL/WGS84 GeoTIFF here:
: https://www.soda-pro.com/help/general-knowledge/linke-turbidity-factor
:: '''''"The Linke turbidity factor (TL, for an air mass equal to 2) is a very convenient approximation to model the atmospheric absorption and scattering of the solar radiation under clear skies. It describes the optical thickness of the atmosphere due to both the absorption by the water vapor and the absorption and scattering by the aerosol particles relative to a dry and clean atmosphere. It summarizes the turbidity of the atmosphere, and hence the attenuation of the direct beam solar radiation. The larger TL, the larger the attenuation of the radiation by the clear sky atmosphere.'''''

:: '''''The TL denotes the transparency of the cloudless atmosphere. If the sky were perfectly dry and clean, TL would be equal to 1. When the sky is deep blue, the TL is just above 1 and still very small. In summer, in Europe, the water vapour is often large and the blue sky is close to white. The TL is larger than 3. In turbid atmosphere, e.g. in polluted cities, the TL is close to 6 - 7."'''''

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get an estimated curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-10-29T05:10:12Z

Hamish: /* Linke turbidity factor values */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. Note that this may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Linke turbidity factor values ==

* See Remund et al. (2003) appendix/table 1 for a list of monthly values for many sites around the world:
: https://www.researchgate.net/publication/42358500_Worldwide_Linke_turbidity_information

* A good description and recipe for calculating it here: (section 3.3.4, page 22)
: https://atmosphere.copernicus.eu/sites/default/files/2022-01/CAMS2_73_2021SC1_D3.2.1_2021_UserGuide_v1.pdf

* Download 5' global dataset for all months of 2010 as LL/WGS84 GeoTIFF here:
: https://www.soda-pro.com/help/general-knowledge/linke-turbidity-factor
:: '''''"The Linke turbidity factor (TL, for an air mass equal to 2) is a very convenient approximation to model the atmospheric absorption and scattering of the solar radiation under clear skies. It describes the optical thickness of the atmosphere due to both the absorption by the water vapor and the absorption and scattering by the aerosol particles relative to a dry and clean atmosphere. It summarizes the turbidity of the atmosphere, and hence the attenuation of the direct beam solar radiation. The larger TL, the larger the attenuation of the radiation by the clear sky atmosphere.

The TL denotes the transparency of the cloudless atmosphere. If the sky were perfectly dry and clean, TL would be equal to 1. When the sky is deep blue, the TL is just above 1 and still very small. In summer, in Europe, the water vapour is often large and the blue sky is close to white. The TL is larger than 3. In turbid atmosphere, e.g. in polluted cities, the TL is close to 6 - 7."'''''

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get an estimated curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-10-29T03:53:26Z

Hamish: /* Automation */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. Note that this may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Linke turbidity factor values ==

* See Remund et al. (2003) appendix/table 1 for a list of monthly values for many sites around the world:
: https://www.researchgate.net/publication/42358500_Worldwide_Linke_turbidity_information

* A good description and recipe for calculating it here: (section 3.3.4, page 22)
: https://atmosphere.copernicus.eu/sites/default/files/2022-01/CAMS2_73_2021SC1_D3.2.1_2021_UserGuide_v1.pdf

* Download 5' global dataset for all months of 2010 as LL/WGS84 GeoTIFF here:
: https://www.soda-pro.com/help/general-knowledge/linke-turbidity-factor

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get an estimated curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T05:08:48Z

Hamish: /* Minimum number of days needed to summarize a solar year */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. Note that this may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get an estimated curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T05:06:37Z

Hamish: /* Tips */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. Note that this may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get a rough curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T05:02:20Z

Hamish: /* Minimum number of days needed to summarize a solar year */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. This may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get a rough curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || ''mid''
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || ''mid''
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || ''mid''
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || ''mid''
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T05:01:58Z

Hamish: /* Optional input maps */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. This may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate ''and'' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get a rough curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || mid
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || mid
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || mid
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || mid
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T04:58:34Z

Hamish: /* Minimum number of days needed to summarize a solar year */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. This may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate 'and' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get a rough curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them.

For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || mid
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || mid
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || mid
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || mid
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T04:57:58Z

Hamish: /* Overview of monthly maps */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. This may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate 'and' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Minimum number of days needed to summarize a solar year ==

You can get a rough curve of daily irradiation throughout the year by running r.sun for e.g. just the middle day of each month, or just for the solstices, equinoxes, and half way between them. For example:

:{| class="wikitable"
|+ Day of year for minimal annual summary
|-
! Day !! DOY !! Comment
|-
| Feb 5th || 036 || mid
|-
| Mar 21st || 080 || equinox
|-
| May 6th || 126 || mid
|-
| Jun 21st || 172 || solstice
|-
| Aug 6th || 218 || mid
|-
| Sept 21st || 264 || equinox
|-
| Nov 6th || 310 || mid
|-
| Dec 21st || 355 || solstice
|}

You may want to duplicate the results from Dec. 21st as day -10 when you fit the curve in your plot.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T04:43:48Z

Hamish: /* Seed maps */ clear up some misinformation

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. This may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Optional input maps ==

Using both input aspect and slope maps (generated from the input DEM) are necessary if you wish to include the effect of reflected light in the global output map.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: potentially buggy''. See {{trac|498}}.

Using horizon seed maps can (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results may not be as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%. Depends on the number of horizon maps. Note for accuracy you might want a large number of pre-calculated horizon maps, but that comes with the downside of it taking a long time and lots of memory to load them off the disk. Calculating horizons in real-time in the CPU may in fact be both more accurate 'and' faster. More tests are needed to clarify this issue.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell. This is only needed if you are running r.sun in a simple XY location. If you are running in a normal session with a known map projection the lat/lon are calculated on the fly.

: The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-29T04:35:24Z

Hamish: /* Tips */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Counterpoint: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.

The <tt>aspin</tt> and <tt>slopein</tt> maps are needed if you wish to include the effect of reflected light (seen in both the <tt>refl_rad</tt> and global radiation output maps). These two maps can be generated from your DEM or DSM with the {{cmd|r.slope.aspect}} module.

The <tt>latin</tt> and <tt>longin</tt> maps are generally only needed if you are running r.sun from a simple XY location. In normal use with a known map projection they are not needed.

If you are in a coastal area you can speed up processing by setting regions of water to NULL in the DEM input raster map. See the {{cmd|r.null}} module. This may reduce the amount of reflected light in the global and <tt>refl_rad</tt> output maps. I suspect that r.sun does not model glint off water at low sun angles, although you may have luck setting water masses to a high albedo with the <tt>albedo</tt> input raster if that is important to you.

== Testing ==

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-11T06:39:36Z

Hamish: /* Overview of monthly maps */ fix dead link

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Dissenting view: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case. The same may be true (or not) for the <tt>latin</tt>, <tt>longin</tt>, <tt>aspin</tt>, and <tt>slopein</tt> pre-calculatated maps.''

=== Testing ===

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [[Time_series#Common_legends_for_many_raster_maps|Common legends for many raster maps]]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-11T04:20:28Z

Hamish: /* Tips */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Dissenting view: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case. The same may be true (or not) for the <tt>latin</tt>, <tt>longin</tt>, <tt>aspin</tt>, and <tt>slopein</tt> pre-calculatated maps.''

=== Testing ===

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [http://grass.osgeo.org/wiki/Time_series#Common_legends_for_many_raster_maps Common legends for many raster maps]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-11T04:17:01Z

Hamish: /* Time step */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Dissenting view: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.''

=== Testing ===

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 '''step=0.5'''

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 '''step=0.25'''

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 '''step=0.05'''

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [http://grass.osgeo.org/wiki/Time_series#Common_legends_for_many_raster_maps Common legends for many raster maps]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-11T04:16:29Z

Hamish: /* Time step */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Dissenting view: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.''

=== Testing ===

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.
: '''''step''' is measured in decimal hours''

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 step=0.5

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 step=0.25

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 step=0.05

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [http://grass.osgeo.org/wiki/Time_series#Common_legends_for_many_raster_maps Common legends for many raster maps]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-11T04:14:05Z

Hamish: /* Tips */

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of {{cmd|r.sun}} may be improved if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.
: ''Dissenting view: modern CPUs are fast enough that the time it takes to read the horizon map off the hard drive may be more than just calculating the horizon on-the-fly, with the added advantage that the on-the-fly horizon is always at exactly the right azimuth. You will need to test to see if r.horizon is better or worse for your own use case.''

=== Testing ===

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 step=0.5

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 step=0.25

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 step=0.05

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [http://grass.osgeo.org/wiki/Time_series#Common_legends_for_many_raster_maps Common legends for many raster maps]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]

R.sun

2025-09-11T02:36:12Z

Hamish: /* Automation */ fix broken link

== Help page ==

* {{cmd|r.sun}} manual page

== Tips ==

The speed of r.sun is much higher if {{cmd|r.horizon}} is used first and the resulting maps are given as input to r.sun. Background: the horizon needs to be computed only one time before, not in every step within r.sun. See the example at the end of the {{cmd|r.sun}} help page.

=== Testing ===

Create an artificial surface containing a Gaussian mound:
{{AddonCmd|r.surf.volcano}} out=gauss method=gaussian kurtosis=1

Overlay some 200m contours to show underlying topography:
r.contour in=gauss out=gauss_200m_contours step=200
d.vect gauss_200m_contours color=white

Set map's color table to highlight detail:
r.colors rad_test.day355.beam col=bcyr -e
d.legend rad_test.day355.beam range=1300,1500

=== Time step ===

The following three images demonstrate the effects of using different time '''step''' parameters.

r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 step=0.5

[[File:Gauss_30min.png|400px|thumb|center|Default 30 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.15minT" day=180 step=0.25

[[File:Gauss_15min.png|400px|thumb|center|15 minute time step over a Gaussian mound.]]

r.sun -s elevin="gauss" glob_rad="rad.global.03minT" day=180 step=0.05

[[File:Gauss_03min.png|400px|thumb|center|3 minute time step over a Gaussian mound.]]

The 3 minute time step takes roughly ten times as long to run as the 30 minute timestep.

== Seed maps ==

Using seed maps can greatly (?, unclear) speed up processing. This is especially important if you will batch process e.g. for every day of the year.

* aspin= and slopein= options: create slope and aspect maps with the {{cmd|r.slope.aspect}} module.
: ''Caution: currently buggy''. See {{trac|498}}.
: Estimated speedup: ?%.

* horizon= and horizonstep= options: Pre-calculate horizon shadows by creating a series of horizon rasters with the {{cmd|r.horizon}} module.
: Results not as smooth as without using this option?? See {{trac|498}}.
: Estimated speedup: ?%.

* latin= and longin= options: Pre-calculate latitudes and longitudes for each raster cell.
The following script will create a raster containing latitude as the raster value:
<source lang="bash">
# create latitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=4 out=elevation.lat
g.remove one
</source>

Creating the longitude map is exactly the same but use column 3 of the {{cmd|m.proj}} output instead of column 4:
<source lang="bash">
# create longitude map (WGS84)
g.region rast=elevation.dem
r.mapcalc one=1
r.out.xyz one | \
cut -f1,2 -d'|' | \
m.proj -oed --quiet | \
sed -e 's/[ \t]/|/g' | \
cut -f1-4 -d'|' | \
r.in.xyz in=- z=3 out=elevation.lon
g.remove one
</source>

: Estimated speedup: 1.3% (in one test).
: See {{trac|498}}.

The above assumes that your projection's ellipsoid is WGS84. If it isn't, use other proj_out= terms with m.proj as required instead of the -o flag.

== Automation ==

It can be tedious to set up and run a long series of r.sun simulations for every day of the year. To automate this we can write a small shell script loop:

<source lang="bash">
for DAY in `seq 1 365` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR
done
</source>

* [https://trac.osgeo.org/grass/browser/grass-addons/grass6/raster/r.sun.tools g.linke_by_day.py] is a small program which will return the Linke coefficient for any day of the year, interpolated from monthly values. You need to edit the script to fill in values appropriate for your study area.

* If you have a multi-core CPU and you'd like to speed things up, here is a small Bourne shell script implementing a poor-man's multi-processing trick:

<source lang="bash">
### r.sun mode 2 loop ###
BEGIN=1
END=365
STEP=1
NUM_CORES=4

for DAY in `seq $BEGIN $STEP $END` ; do
DAY_STR=`echo $DAY | awk '{printf("%.03d", $1)}'`
echo "Processing day $DAY_STR at `date` ..."

LINKE="`g.linke_by_day.py $DAY`"

CMD="r.sun -s elevin=elevation.dem day=$DAY lin=$LINKE step=0.05 \
beam_rad=rad_beam.$DAY_STR diff_rad=rad_diffuse.$DAY_STR \
refl_rad=rad_reflected.$DAY_STR glob_rad=rad_global.$DAY_STR \
insol_time=rad_insol_time.$DAY_STR --quiet"

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] || [ "$DAY" = "$END" ] ; then
# stall to let the background jobs finish
$CMD
sleep 2
wait
#while [ `pgrep -c r.sun` -ne 0 ] ; do
# sleep 5
#done
else
$CMD &
fi
done
wait # wait for background jobs to finish to avoid race conditions
</source>

* For a r.sun Mode 1 loop, see the [[Parallelizing Scripts]] page.

== Overview of monthly maps ==

Given there are mean monthly global irradiation maps calculated, an overview map of all 12 maps can be drawn along with their legends (see also [http://grass.osgeo.org/wiki/Time_series#Common_legends_for_many_raster_maps Common legends for many raster maps]). The script assumes a common naming convention for all 12 monthly average maps (i.e. global_rad_avg.jan, global_rad_avg.feb, global_rad_avg.mar, etc.).

<source lang="bash">
#!/bash/sh
# script to draw monthly (mean) global solar irradiation maps in a 3x4 matrix

# set wide aspect ratio (16:9, e.g. 1366 width x 768 height)
d.monsize setm=x0 setw=1366 seth=768

# split in 12 frames
d.split.frame frames=12

# preferred font?
d.font DroidSans

for FRAME in "uno dec December" "dos jan January" "tres feb February" \
"cuatro mar March" "cinco apr April" "seis may May" \
"siete jun June" "ocho jul July" "nueve aug August" \
"diez sep September" "once oct October" "doce nov November"

do

# parse "${FRAMES_STR}" and set positional parameters
set -- $FRAME ; echo $1 $2 $3

# select FRAME
d.frame -s "${1}"

# draw the map
d.rast global_rad_avg."${2}"

# draw label on the left and vertically
d.text text="${3}" size=10 color=50:50:50 at=9,25 rotation=90

# draw legend
d.legend global_rad_avg."${2}" -s at=10,90,89,92

done
</source>

[[File:Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece.png|400px|thumb|center|Mean monthly global irradiation maps (n=38.918, s=38.531, w=21.793, e=22.436)]]
Same image also [http://grass.osgeo.org/grass-wiki/images/Mean_monthly_global_rad_maps_over_a_study_area_in_central_Greece_700px_height.jpg here (a bit of higher resolution)]

== Development ==

* Ongoing trials documented in [https://trac.osgeo.org/grass/ticket/498 trac ticket #498]

=== OpenCL ===
As part of the Google Summer of Code 2010 an {{wikipedia|OpenCL}} version has been written that allows ''r.sun'' to run on {{wikipedia|GPU}}s. This provides a massive speedup in processing time. The merging of this code into GRASS 7's r.sun is forthcoming (HB).
* To get your hands on the code now see http://github.com/mailseth/OpenCL-integration-for-GRASS---GDAL (functional prototype now available)
* To use OpenCL GPU functionality you'll want a graphics card like the ATI HD5770 or nVidia GTX260 or GTX460, or newer.
* OpenCL also allows for running on multi-core (or multi-processor) CPUs, for systems without GPUs. You can run OpenCL via multicore *or* GPU, but not both on the same job at the same time.
* Seth wrote: ''The OpenCL version of r.sun runs over 20x faster than the original version on my machine (2.26 GHz Mac Pro vs. GeForce GTX 285). However, it is hampered by the low memory on your GPU, so you may need to partition your raster.''
* You can now ./configure GRASS 7 with OpenCL support using:
--with-opencl

== ToDo ==

* Merge in OpenCL code
* Add support for [[OpenMP]] parallelization (if OpenCL doesn't make that redundant)

== References ==

Usage of r.sun in the literature (selected references):

* Agugiaro, G., Remondino, F., Stevanato, G., De Filippi, R., Furlanello, C., 2011. Estimation of solar radiation on building roofs in mountainous areas. International Archives of Photogrammetry 38. PDF: http://www.pf.bv.tum.de/isprs/pia11/pub/PIA11_Agugiaro_et_al.pdf
* H.T. Nguyen and J.M. Pearce, "[http://dx.doi.org/10.1016/j.solener.2010.02.009 Estimating Potential Photovoltaic Yield with r.sun and the Open Source Geographical Resources Analysis Support System]" Solar Energy 84, pp. 831-843, 2010. [http://mtu.academia.edu/JoshuaPearce/Papers/1540702/Estimating_Potential_Photovoltaic_Yield_with_r.sun_and_the_Open_Source_Geographical_Resources_Analysis_Support_System Open Access] (see also [http://www.appropedia.org/An_open_source_simulation_of_photovoltaic_yield_with_r.sun_over_large_regions Link to related external wiki page])
* Hofierka, J., Kanuk, J., 2009. Assessment of photovoltaic potential in urban areas using open-source solar radiation tools. Renewable Energy 34, 2206-2214.
* Huld, T.A., Šúri, M., Dunlop, E.D., Micale, F., 2006. Estimating average daytime and daily temperature profiles within Europe. Environmental Modelling & Software 21, 1650-1661.
* Šúri, M., Hofierka, J., 2004. A new GIS‐based solar radiation model and its application to photovoltaic assessments. Transactions in GIS 8, 175-190.

[[Category: Documentation]]
[[Category: Parallelization]]