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

GRASS GIS GRASS GIS developer
1

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

(attachments)

rsun_no_pj.diff (2.4 KB)

2

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

3

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

4

Seth wrote:

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.

Jaro's patch gives the same exact result as mine, but his is faster,
simpler, and more in sync with the historic code. I still have to study
Thomas's method a bit more to fully understand, but for now my observer's
money would be on taking Jaro's patch.

cheers,
Hamish

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);

and yet mine does as well, but does not show the rotation.. ?

> 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:
>> 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?