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

GDAL · June 21, 2013, 3:17pm

#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: | Platform: Linux
Cpu: x86-64 |
---------------------------+------------------------------------------------
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

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

GDAL · June 22, 2013, 4:09am

#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 |
---------------------------+------------------------------------------------
Changes (by hamish):

* keywords: => r.sun

Comment:

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

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

GDAL · September 9, 2013, 3:20am

#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 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

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

GDAL · September 9, 2013, 10:14am

#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 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

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

GDAL · September 9, 2013, 10:18am

#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):

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)

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

GDAL · September 15, 2013, 3:36am

#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, 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

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

GDAL · September 15, 2013, 3:54am

#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):

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

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

GDAL · September 15, 2013, 11:32pm

#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 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

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

GDAL · September 16, 2013, 2:40am

#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):

Replying to [comment:6 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.

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

GDAL · September 25, 2013, 5:13am

#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>

GDAL · September 25, 2013, 5:49am

#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):

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

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