[GRASS-dev] r.sun and pj_do_proj() calls

[Hamish]

Hi,

As suggested by Jaro, by
- commenting out all the pj_do_proj() stuff, and
- re-enabling sin(sunVarGeom->sunAzimuthAngle), and
- using the -s shading flag, and
- not using horizon maps,

I get the same exact output values as with my patch (attached), which
replaces the proper pj_do_proj() reprojection with a simple cosine
pseudo-projection.

Both Jaro's and my patches give slightly different results than the
current SVN version based on pj_do_proj(). Looking at the test results
from using a Gaussian mound the pj_() version seems to be rotated counter-
clockwise by a little bit compared to the others. My best guess is that
this is due to the 0.9 degree difference between grid-north and
true-north at this site (`g.region -n`), as it is the pj_do_proj()
version which seems to be slightly askew, the others seem to be 
symmetric. But perhaps that rotation is more than 1 deg?


# spearfish dataset
g.region -dp
r.surf.volcano out=gauss method=gaussian kurtosis=1


# current r.sun in grass 6.5svn ( using pj_do_proj() )
time r.sun -s elevin="gauss" glob_rad="rad.global.30minT.65svn" day=180
step=0.5 real    3m19.480s
  user    3m6.108s
  sys     0m2.500s


# Jaro's patch ( revert back to setting from sunAzimuthAngle )
time r.sun -s elevin="gauss" glob_rad="rad.global.30minT.sunAzimuthAngle"
day=180 step=0.5 real    2m58.834s
  user    2m50.555s
  sys     0m1.188s


# My patch ( attached; replace pj_do_proj() with lon*cos(lat) scaling )
time r.sun -s elevin="gauss" glob_rad="rad.global.30minT.coslat" day=180
step=0.5 real    3m12.652s
  user    2m52.727s
  sys     0m1.424s



# compare results
for map in 65svn sunAzimuthAngle coslat ; do
  echo "[$map]"
  r.univar rad.global.30minT.$map -g | grep 'mean=\|stddev=\|sum='
  echo
done

[65svn]
mean=8788.2168789737
stddev=40.3700050839272
sum=2657714972.10546875

[sunAzimuthAngle]
mean=8788.04781049377
stddev=40.3848632592666
sum=2657663842.75390625

[coslat]
mean=8788.04781049377
stddev=40.3848632592666
sum=2657663842.75390625


# view
d.erase
for map in 65svn sunAzimuthAngle coslat ; do
   echo "[$map]"; d.rast rad.global.30minT.$map
   d.title -s rad.global.30minT.$map | d.text
   read
done


(the Gaussian mound also shows why a small step size like 0.05 is so
important, which may finally be made practical by GPU acceleration. see
  http://grass.osgeo.org/wiki/r.sun#Time_step )


regards,
Hamish
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rsun_no_pj.diff
Type: text/x-diff
Size: 2457 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-dev/attachments/20100801/7191b2a8/rsun_no_pj-0001.bin