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
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
