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

[Jaro Hofierka]

Hi Hamish,

Many thanks for your great work!
I expected some kind of skewness in the svn version, because it always 
seeks a point in one direction:

>     delt_lat = -0.0001 * cos(inputAngle);  /* Arbitrary small distance in latitude */
>     delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);

My idea is that perhaps the easiest way to deal with the meridian 
convergence would be to do all calculations in lat-lon, even for 
cartesian systems.
Even now we call pj_do_proj() to get lat-lon for every grid point in 
main.c, so it might be more convenient to get the lat-lon values in the 
beginning and then do everything in lat-lon.

Kind regards,

Jaro


Hamish wrote:
> 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