[GRASS-dev] [GRASS GIS] #2014: r.sun using EPSG:3031 projection gives strange results

[GRASS GIS]

#2014: r.sun using EPSG:3031 projection gives strange results
---------------------------+------------------------------------------------
 Reporter:  pierreroudier  |       Owner:  grass-dev@…              
     Type:  defect         |      Status:  new                      
 Priority:  normal         |   Milestone:  7.0.0                    
Component:  Raster         |     Version:  svn-trunk                
 Keywords:  r.sun          |    Platform:  Linux                    
      Cpu:  x86-64         |  
---------------------------+------------------------------------------------

Comment(by hamish):

 Hi, fwiw after fixing a small bug in r.sunmask I was able to write a
 script to show the sun's position through the day on a polar graph.

 Plot is for Jan 1st 00:00-23:45 (gap in the plot is 23:45-24:00) around
 Ross Island in Antarctica (McMurdo & Scott Bases) in a polar stereographic
 location (epsg:3286). Red "+" are the 15 minute marks. So the more sun in
 the south effect seems to be a combination of midnight shade on the north
 slope, sun peeking over the ridges onto the south slope at midday, and the
 cumulative time the Sun spends in the south at this time of year- the
 azimuth is not strictly going around at 15deg/hour.

 d.histogram of a r.slope.aspect slope map shows most of Erebus between 4
 and 17 deg, very little over 25 deg, and a max of 35deg. So midday sun at
 35.5 deg above the horizon would generally get to the south slope, even if
 the shallow angle meant that little of it would be hitting it with much
 impact.

 so it's all about the cumulative timing & the tortoise beats the hare.


 Hamish


 -----
 The r.sunmask module can calculate the position of the Sun in the sky
 throughout the day:
 {{{
 for HOUR in `seq -f'%02.0f' 0 23` ; do
   for MINUTE in 00 15 30 45 ; do
     r.sunmask -sg elev=ramp200dem_wgs_v2.ross_island out=dummy \
       year=2013 month=1 day=1 hour=$HOUR minute=$MINUTE second=0 \
       timezone=12 \
      > "solar_pos.hour_$HOUR.$MINUTE.txt"
   done
 done

 grep sun solar_pos.hour_* | grep -v set | grep -v rise > solar_pos.day
 grep above solar_pos.day | cut -f2 -d= > angleabove.txt
 grep azim solar_pos.day | cut -f2 -d= > azimuth.txt

 % matlab or octave
 load azimuth.txt
 load angleabove.txt
 polar(azimuth * pi/180, angleabove)
 hold on
 polar(azimuth * pi/180, angleabove, 'r+')
 set(gca, 'View', [90 -90])  % rotate so 0 deg as north-up not +x axis
 }}}

-- 
Ticket URL: <https://trac.osgeo.org/grass/ticket/2014#comment:9>
GRASS GIS <http://grass.osgeo.org>