Should I go with Hamish's patch in the previous email, or is there some ongoing work here? I think I'm getting ready to push to finish OpenCL r.sun.
~Seth

On Aug 1, 2010, at 12:07 AM, Jaro Hofierka wrote:

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
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-dev

Reply via email to