[toc]

Recently, in order to nicely plan ahead for a birthday lunch at the agriturismo Malga Brigolina (a farm-restaurant near Sopramonte di Trento, Italy, at 1,000m a.s.l. in the Southern Alps), friends of mine asked me the day before:

“Will the place be sunny at lunch time for making a nice walk?“

[well, the weather was close to clear sky conditions but mountains are high here and cast long shadows in winter time].

A rather easy task I thought, so I got my tools ready since that was an occasion to verify the predictions with some photos! Thanks to the new EU-DEM at 25m I was able to perform the computations right away in a metric system rather than dealing with degree in LatLong.

Direct sunlight can be assessed from the beam radiation map of GRASS GIS’ r.sun when running it for a specific day and time. But even easier, there is a new Addon for GRASS GIS 7 which calculates right away time series of insolation maps given start/stop timestamps and a time step: r.sun.hourly. This shrinks the overall effort to almost nothing.

Creating the direct sunlight maps

The first step is to calculate where direct sunlight reaches the ground. Here the input elevation map is the European “eu_dem_25”, while the output is the beam radiation for a certain day (15 Dec 2013 is DOY 349).

Important hint: the computational region must be large enough to east/south/west to capure the cast shadow effects of all relevant surrounding mountains.

I let calculations start at 8am and finish at 5pm which an hourly time step. The authors have kindly parallelized r.sun.hourly, so I let it run on four processors simultaneously to speed up:

# calculate DOY (day-of-year) from given date:
date -d 2013-12-15 +%j
349

# calculate beam radiation maps for a given time period
# (note: minutes are to be given in decimal time: 30min = 0.5)
r.sun.hourly elev_in=eu_dem_25 beam_rad=trento_beam_doy349 \
      start_time=8 end_time=17 time_step=1 day=349 year=2013 \
      nprocs=4

The ten resulting maps contain the beam radiation for each pixel considering the cast shadow effects of the pixel-surrounding mountains. However, the question was not to calculate irradiance raster maps in Wh/m^2 but simply “sun-yes” or “sun-no”. So a subsequent filtering had to be applied. Each beam radiation map was filtered: if pixel value equal to 0 then “sun-no”, otherwise “sun-yes” (what my friends wanted to achieve; effectively a conversion into a binary map).
Best done in a simple shell script loop:

[Edit 30 Dec 2013: thanks to Anna you can simplify below loop to the r.colors call thanks to the new -b flag for binary output in r.sun.hourly!]

for map in `g.mlist rast pattern="trento_beam_doy349*"` ; do
    # rename current map to tmp
    g.rename rast=$map,$map.tmp
    # filter and save with original name
    r.mapcalc "$map = if($map.tmp == 0, null(), 1)"
    # colorize the binary map
    echo "1 yellow" | r.colors $map rules=-
    # remove cruft
    g.remove rast=$map.tmp
done

As a result we got ten binary maps, ideal for using them as overlay with shaded DEMs or OpenStreetMap layers. The areas exposed to direct sunlight are shown in yellow.

Time to look at some details: So, is Malga Brigolina in sunlight at lunch time?

Situation at 12pm (noon): predicted that the restaurant is still in shadow – confirmed in the photo:

Situation at 1:30/2:00pm: sun is getting closer to the Malga, as confirmed in photo (note that the left photo is 20min ahead of the map). The small street in the right photo is still in the shadow as predicted in the map):

(click to enlarge)

Situation at 3:00pm: sun here and there at Malga Brigolina:

Sunlight map blended with OpenStreetmap layer (r.blend + r.composite) and draped over DEM in wxNVIZ of GRASS GIS 7 (click to enlarge). The sunday walk path around Malga Brigolina is the blue/red vector line shown in the view center.

Situation at 4:00pm: we are close to sunset in Trentino… view towards the Rotaliana (Mezzocorona, S. Michele all’Adige), last sunlit summits also seen in photo:

(click to enlarge)

Outcome

The resulting sunlight/shadow map appear to match nicely realty. Perhaps r.sun.hourly should get an additional flag to generate the binary “sun-yes” – “sun-no” maps directly.

Direct sunlight zones (yellow, 15 Dec 2013, 3pm): Trento with Monte Bondone, Paganella, Marzolan, Lago di Caldonazzo, Lago di Levico and surroundings (click to enlarge)

GRASS GIS usage note

The wxGUI settings were as simple as this (note the transparency values for the various layers):

Data sources:

• Shaded maps produced using Copernicus data and information funded by the European Union – EU-DEM layers
• © OpenStreetMap contributors