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

[Seth Price]

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