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
Index: raster/r.sun/rsunlib.c
===================================================================
--- raster/r.sun/rsunlib.c (revision 42938)
+++ raster/r.sun/rsunlib.c (working copy)
@@ -188,14 +188,11 @@
struct GridGeometry *gridGeom, double latitude, double longitude)
{
double pom, xpom, ypom;
-
double costimeAngle;
double lum_Lx, lum_Ly;
-
- double newLatitude, newLongitude;
double inputAngle;
double delt_lat, delt_lon;
- double delt_east, delt_nor;
+ double delt_lat_m, delt_lon_m;
double delt_dist;
@@ -247,7 +244,6 @@
}
-
if (sunVarGeom->solarAzimuth < 0.5 * M_PI)
sunVarGeom->sunAzimuthAngle = 0.5 * M_PI - sunVarGeom->solarAzimuth;
else
@@ -257,33 +253,22 @@
inputAngle = sunVarGeom->sunAzimuthAngle + pihalf;
inputAngle = (inputAngle >= pi2) ? inputAngle - pi2 : inputAngle;
-
+ /* 1852m * 60 * 0.0001rad * 180/pi= 636.67m */
delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
- newLatitude = (latitude + delt_lat) * rad2deg;
- newLongitude = (longitude + delt_lon) * rad2deg;
+ delt_lat_m = delt_lat * (180/M_PI) * 1852*60;
+ delt_lon_m = delt_lon * (180/M_PI) * 1852*60 * cos(latitude);
+ delt_dist = sqrt(delt_lat_m * delt_lat_m + delt_lon_m * delt_lon_m);
+/*
+ sunVarGeom->stepsinangle = gridGeom->stepxy * sin(sunVarGeom->sunAzimuthAngle);
+ sunVarGeom->stepcosangle = gridGeom->stepxy * cos(sunVarGeom->sunAzimuthAngle);
+*/
- if ((G_projection() != PROJECTION_LL)) {
- if (pj_do_proj(&newLongitude, &newLatitude, &oproj, &iproj) < 0) {
- G_fatal_error("Error in pj_do_proj");
- }
- }
+ sunVarGeom->stepsinangle = gridGeom->stepxy * delt_lat_m / delt_dist;
+ sunVarGeom->stepcosangle = gridGeom->stepxy * delt_lon_m / delt_dist;
- delt_east = newLongitude - gridGeom->xp;
- delt_nor = newLatitude - gridGeom->yp;
-
- delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor);
-
-
- sunVarGeom->stepsinangle = gridGeom->stepxy * delt_nor / delt_dist;
- sunVarGeom->stepcosangle = gridGeom->stepxy * delt_east / delt_dist;
-
- /*
- sunVarGeom->stepsinangle = stepxy * sin(sunVarGeom->sunAzimuthAngle);
- sunVarGeom->stepcosangle = stepxy * cos(sunVarGeom->sunAzimuthAngle);
- */
sunVarGeom->tanSolarAltitude = tan(sunVarGeom->solarAltitude);
return;
_______________________________________________
grass-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-dev
