r.sun using EPSG:3031 projection gives strange results

[pierreroudier]

I am working on Antarctica data, projected in Antarctic Polar Stereographic (EPSG:3031 [0][1]). This projection puts the South Pole in the "center" of the map.

I have strange results in the Ross Sea Region using r.sun from GRASS 7 compiled from trunk: According to the r.sun results, the south facing slope receive more radiation than the north facing one, which doesn't add up.

[0] ​http://spatialreference.org/ref/epsg/3031/

[1] ​http://nsidc.org/data/atlas/epsg_3031.html

[hamish]

Hi,

I saw the post on the mailing list and will have a look at it -- it might take me a little while to get to it though. You might try to compare with the version of r.sun in GRASS 6.4.3svn, as in GRASS 7 some map projection calculations have been moved out of the main loop to allow the GPU/multiprocessor acceleration to work.

It won't necessarily fix it, but it may be useful to test.

A couple other things to look at: is the sun over or under the horizon when it happens (there may be no direct sun but still some diffuse radiation?), and if PROJ.4 is handling the projection properly (projections not "north up" (which way is north in a polar stereographic?), non-Greenwich prime meridians etc.)

in general, polar stereographic is known to work, at least with v.proj:

​http://grassold.osgeo.org/screenshots/images/day_on_earth_S.png

Hamish

[pierreroudier]

test on Ross Sea Region

[pierreroudier]

Hi Hamish,

Sorry it took me so long to get back to you, but finally got some time for further testing,

If you look at the attached file, you can see that it seems that r.sun is integrating the latitude and altitude effects OK, but the shadowing seem to be the wrong way round: on the Ross Sea Region pictured, the South Pole is on the "top" of the map, so faces looking "up" should receive less solar energy than faces looking "down".

Does that make sense?

Cheers,

Pierre

[hamish]

Hi Pierre,

could you try to make an artificial Gaussian surface using r.surf.volcano as described in ticket #498 and/or ​http://grasswiki.osgeo.org/wiki/R.sun#Testing ?

it helps make the issues obvious.

do you know of a readily available DEM of the area we could test with? (~169E, 73.5S ?)

draping the sunlight over the DEM in NVIZ can help.

direct rad only? what's the exact r.sun command line used?

what happens if you try another custom polar stereographic with +lon_0=170 as "up"? and then +lon_0 as 90E or 90W for "up"? or a custom local tmerc centred on your study site?

regards, Hamish

[hamish]

btw, was the screenshot made using GRASS? (if so I'd like to do something about fixing whatever's responsible for that .000000 %f formatted output seen in the legend)

[hamish]

Hi, I had a try using Mt. Erebus as the shadowing body & r.sun in mode 2 (integrated over the entire day),

north:      -1280000
south:      -1372500
west:       240000
east:       339500
nsres:      100
ewres:      100

and cubic-reprojected etopo1 for elevation. (can probably dig out some lidar for Ross Island from someone in our dept; 100m DEM from Colorado is currently offline due to flooding)

test setup for grass 6.5 looked like this:

DAYS="001 030 045 060 070 080"
for DAY in $DAYS ; do
   r.sun -s elevin=etopo1.ross_island step=0.05 \
      beam_rad=erebus.direct.$DAY diff=erebus.diffuse.$DAY day=$DAY &
done
wait

DAYS="085 090 095 100 105 110 115"
for DAY in $DAYS ; do
   r.sun -s elevin=etopo1.ross_island step=0.05 \
      beam_rad=erebus.direct.$DAY diff=erebus.diffuse.$DAY day=$DAY &
done

for map in `g.mlist pat=erebus.d*` ; do
   r.colors $map col=bcyr
done

day 115 is near last-light of the autumn, not exactly the same sunrise/set as these values for McMurdo, but pretty close & I'm not sure which definition they're using.

​http://www.timeanddate.com/worldclock/astronomy.html?n=1032&month=4&year=2013&obj=sun&afl=-11&day=1

There's an obvious low-sun-angle bug in the diffuse output map, which I'll file in a separate report, and obviously the default albedo and linke values in the module won't be much good, but the direct beam results seem ok if I look at the individual days. I'm still waiting for days 1 & 30 to finish, but looking at the result for day 45 makes it seem like for summertime there might be some counter-intuitive irradiation pattern, with more light on the true-south side of the mountain. If day 1 or 30 show the same, the next thing would be to run r.sun in mode 1 quarter-hourly over that one day and see if that can make more sense.

todo: test with grass7. The main differences in the versions is that by default shadowing effect of the DEM is turned on in G7, while in grass 6.x you needed to use the -s flag; and the lat/lon reprojection code is backed out of the main loop for speed and to make parallelization/GPU assist easier.

Hamish

[hamish]

yeah, looking at the results for day=1 and 30 there's significantly more light showing on the true-south side of Erebus at that time of year. To test if that's a counter-intuitive cumulative time-spent effect or not, next things to try I think are to plot the sun's position through the day, and try it in another Antarctic stereographic projection that puts grid-north and true-north in alignment and see if we get the same result.

Hamish

[pierreroudier]

Hamish,

The original command I used was:

## annual beam radiance
DEM=rsr_dem

# setup region, and re-compute slope in deg.
 g.region rast=$DEM res=200 -ap
#g.region res=200 n=-567200 s=-634800 w=-3400 e=117600 -ap

r.slope.aspect --o elev=$DEM slope=slope_deg aspect=aspect

# setup: no accomodation for Linke, using defaults
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` ..."
echo "Processing day $DAY_STR at `date` ...\n" >> r.sun_progress.txt

CMD="r.sun --o --quiet \
	 elev_in=$DEM slope_in=slope_deg asp_in=aspect \
         step=1.0 \
	 day=$DAY \
         beam_rad=beam_rad.$DAY_STR \
         insol_time=insol_time.$DAY_STR "

# poor man's multi-threading for a multi-core CPU
MODULUS=`echo "$DAY $STEP $NUM_CORES" | awk '{print $1 % ($2 * $3)}'`
if [ "$MODULUS" = "$STEP" ] ; then
   # stall to let the background jobs finish
   $CMD
   sleep 2
else
   $CMD &
fi
done


# combine into annual sum
r.series --o in=`g.mlist type=rast pattern="beam_rad.*" sep=","` out=beam_rad_annual method=sum 
r.series --o in=`g.mlist type=rast pattern="insol_time.*" sep=","` out=insol_time_annual method=sum

# convert to MJ/sq. m
# r.sun mode 2 * (3600 / 100000) ==> MJ/sq. m
r.mapcalc --o "beam_rad_annual_mj = beam_rad_annual * (3600.0 / 100000.0)"
r.colors.stddev beam_rad_annual_mj


# clean-up
# g.mremove -f rast=beam_rad.*
# g.mremove -f rast=insol_time.*

Happy to assist with the testing on grass7. Which Antarctic projection would you suggests? I got no projection in mind that would suit our needs.

Pierre

BTW the screenshot was generated in QGIS, not GRASS, so that's one less bug to fix :)

[hamish]

Replying to hamish:

yeah, looking at the results for day=1 and 30 there's significantly more light showing on the true-south side of Erebus at that time of year. To test if that's a counter-intuitive cumulative time-spent effect or not, next things to try I think are to plot the sun's position through the day, and try it in another Antarctic stereographic projection that puts grid-north and true-north in alignment and see if we get the same result.

Just tried with epsg:3286 (ant. polar stereo +lon_0=165) and the results are the same. So I think what we're seeing is that at in mid-summer the low sun angle from the south at midnight is shadowing the north side of Erebus, while at midday the sun is high enough in the sky to get a bit of incident light onto the southern slope of the mountain. Net result: in mid-summer there seems to be more sunlight reaching the southern side of the mountains in Antarctica.

Near first/last light in spring/autumn with the sun just peeking over the northern horizon the shadow effects are as you'd expect.

It would be good to try with more realistic values for albedo (since snow) and Linke turbidity (since the beams are passing through a lot of atmosphere at these angles), and also to do the fine resolution sum for January 1st by hour or finer to confirm that the tortoise beats the hare.

Day 045 seems to be where the crossover happens from more sun reaching the south to more reaching the north side.

regards, Hamish

ps- glad to see that 'g.region -n' is working in both projections, but d.grid is only working in epsg:3286 so had to use v.mkgrid + v.proj to get the lat/lon (1 deg x 20 minutes) overlay grid there.

[hamish]

Same as previous with 15minute marks

[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

[hamish]

(radius of the plot line is angle of the Sun above the horizon, theta is compass angle, north is 0 deg, east is 90, ..)

[martinl]

What is state of this ticket?

[martinl]

What is status of this issue?