[GRASS-user] Question about r.sun2 and r.horizon

GRASS GRASS user
1

Dear GRASS-users,
Dear Mr. Huld,

Please accept my apologies for this long post and for double-posting
this. My "felis" account in the university does not seem to deliver this
messages in the "jrc.it" mail-server.

I have following question about r.sun2 (and r.horizon):

My colleague Dimitris pointed-out that the global_radiation maps
produced with r.sun2 show a higher radiation (for the same elevation
more or less) in north expositions comparing with southern expositions.
I confirm this behaviour after careful examination of the maps I
produced (if interested I could send you a screen-shot or more with
fragment(s) from
our average global_rad map in kJ).

We expect the exact opposite effect on the globa_rad maps. What could be
my mistake?

###

Details about the map production using the following data:

*srtm3

*slope & aspect produced with r.slope.aspect

*linke turbidity = 5.2
(average value for the region

N: 39:05:38.098741N S: 38:23:35.098741N Res: 0:00:03
E: 22:26:13.344612E W: 21:21:01.344612E Res: 0:00:03

and for the period April to October = julian days 92 to 304)

*albedo = 0.23

I produced global_radiation maps using r.sun2 based on the following

1.

r.horizon elevin=srtm3_stel step=45 horizon=horizon_45 coord=4275945

*** Question about r.horizon
************************************************************************************

The above command gives horizon maps like hortest.45_0, hortest.45_1,
etc. Following the r.sun2 expects something like hortest.45_000,
hortest.45_001, etc.

If r.sun2 expects a specific prefix-pattern for r.horizon maps, why does
r.horizon name its ouput maps different than what r.sun2 expects?

**************************************************************************************************************************

2.
(using r.sun2 in a loop)

for i in `seq 92 304`; do r.sun2 -s --o elevin=srtm3_stel
aspin=aspect_stel slopein=slope_stel alb=0.23 lin=5.2 lat=4275945
horizon=horizon.45 horizonstep=45 glob_rad=global_rad_$i day=$i; done

3.

for i in global_rad_*; do r.mapcalc global_daily_avg=global_daily_avg +
$i; echo "...adding map $i"; done
        
4.

r.mapcalc global_daily_avg=global_daily_avg/212

5.

r.mapcalc global_daily_avg_kJ="int(global_daily_avg * 3.6)"

Kind regards

Nikos

2

On Sun, 2008-06-15 at 22:30 +0200, Nikos Alexandris wrote:

1.

r.horizon elevin=srtm3_stel step=45 horizon=horizon_45 coord=4275945

# This is a bug and will be fixed.

2.
(using r.sun2 in a loop)

for i in `seq 92 304`; do r.sun2 -s --o elevin=srtm3_stel
aspin=aspect_stel slopein=slope_stel alb=0.23 lin=5.2 lat=4275945
horizon=horizon.45 horizonstep=45 glob_rad=global_rad_$i day=$i; done

I created latitude and longitude maps, re-run r.sun2 and compared a day
(on derived from the above command and another one derived from the
command below). The results seem rational now.

r.sun2 -s elevin=srtm3_for_global_rad aspin=aspect_for_global_rad
slopein=slope_for_global_rad alb=0.23 lin=5.2
latin=latitude_for_global_rad longin=longitude_for_global_rad
horizon=horizon.45 horizonstep=45 glob_rad=global_rad_testday92 day=92

So the question is now, why when using a fixed lat value, the northern
expositions have a higher radiation than the southern expositions?

3

> r.horizon elevin=srtm3_stel step=45 horizon=horizon_45 coord=4275945

coord= should be x,y not just a single number.

So the parser knows to expect two numbers, in the option setup code should include the line:
  parm.coord->key_desc = "east,north";

> for i in `seq 92 304`; do
> r.sun2 -s --o elevin=srtm3_stel aspin=aspect_stel \
> slopein=slope_stel alb=0.23 lin=5.2 lat=4275945 \
> horizon=horizon.45 horizonstep=45 \
> glob_rad=global_rad_$i day=$i; done

lat= must be in the range 0-90. it will be converted to radians and so does not take projected northings as input.

the parser option to ensure the 0-90 range is:
  parm.lat->options = "0-90";

although it might be nicer to allow DDD:MM:SS.SSS strings (e.g. "23:55N") then use G_lat_scan() to convert to a float. check the return value of that function to see if it is out of range and a G_fatal_error() is needed.

but as far as I can tell looking through the r.sun2 code, the *lt string is only ever used when writing out the map's metadata/history. (??)

Hamish

4

On Mon, 2008-06-16 at 07:03 -0700, Hamish wrote:

> > r.horizon elevin=srtm3_stel step=45 horizon=horizon_45 coord=4275945

coord= should be x,y not just a single number.

So the parser knows to expect two numbers, in the option setup code should include the line:
  parm.coord->key_desc = "east,north";

I remember using the centre x,y in the beginning and changing it later
only to x because I read the r.sun2 help which is:

[...]
lat A single value of latitude
[...]

> > for i in `seq 92 304`; do
> > r.sun2 -s --o elevin=srtm3_stel aspin=aspect_stel \
> > slopein=slope_stel alb=0.23 lin=5.2 lat=4275945 \
> > horizon=horizon.45 horizonstep=45 \
> > glob_rad=global_rad_$i day=$i; done

lat= must be in the range 0-90. it will be converted to radians and so does not take projected northings as input.

You give me more homework :wink: I didn't really checked that and the loop
is running the last 2 hours!!!

the parser option to ensure the 0-90 range is:
  parm.lat->options = "0-90";

although it might be nicer to allow DDD:MM:SS.SSS strings (e.g. "23:55N") then use G_lat_scan() to convert to a float. check the return value of that function to see if it is out of range and a G_fatal_error() is needed.

but as far as I can tell looking through the r.sun2 code, the *lt string is only ever used when writing out the map's metadata/history. (??)

Hamish

5

Nikos:

r.horizon elevin=srtm3_stel step=45 horizon=horizon_45 coord=4275945

Hamish:

> coord= should be x,y not just a single number.
>
> So the parser knows to expect two numbers, in the
> option setup code should include the line:
> parm.coord->key_desc = "east,north";

(added in SVN)

I remember using the centre x,y in the beginning and
changing it later only to x because I read the r.sun2 help which
is:

[...]
lat A single value of latitude
[...]

the r.horizon code expects %lf,%lf. see main.c line 343.
http://trac.osgeo.org/grass/browser/grass-addons/raster/r.sun_horizon/r.horizon/main.c#L343

> > for i in `seq 92 304`; do
> > r.sun2 -s --o elevin=srtm3_stel aspin=aspect_stel \
> > slopein=slope_stel alb=0.23 lin=5.2 lat=4275945 \
> > horizon=horizon.45 horizonstep=45 \
> > glob_rad=global_rad_$i day=$i; done

Hamish:

> lat= must be in the range 0-90. it will be converted
> to radians and so does not take projected northings as
> input.

....

> but as far as I can tell looking through the r.sun2
> code, the *lt string is only ever used when writing out the
> map's metadata/history. (??)

I added a comment about G_scan_lat() into the code, but the lat= answer is still unused by the code AFAICT. I am not sure of the authors' intentions so left it as-is.

don't feel bad, it is very important to find these things..

Hamish

6

Just reporting how I finally calculated the global_rad map.

r.slope.aspect elev=srtm3_for_global_rad slope=slope_for_global_rad

aspect=aspect_for_global_rad format=degrees

i.latitude input=srtm3_for_global_rad latitude=latitude_for_global_rad

i.longitude input=srtm3_for_global_rad

longitude=longitude_for_global_rad

r.horizon elevin=srtm3_for_global_rad step=45 horizon=horizon45 ###

This gives radians as an output!

for i in `seq 92 304`; do r.sun2 -s --o elevin=srtm3_for_global_rad

aspin=aspect_for_global_rad slopein=slope_for_global_rad alb=0.23
lin=5.2 horizon=horizon45 horizonstep=45 glob_rad=global_rad_$i day=
$i; done

### After Hamish' comments and a bit of testing I found no difference
when using different values on the "lat=" parameter. So it really stays
out of the game (?)

r.mapcalc global_rad_sum=0 && for i in `seq 92 304`; do r.mapcalc

global_rad_sum=global_rad_sum + global_rad_$i; done

r.mapcalc global_daily_avg=global_rad_sum / 212

r.mapcalc global_daily_avg_kJ="int(global_daily_avg * 3.6)"

r.mapcalc formula="global_daily_avg_kJ / 2450 * (0.025000 *

temperature_int + 0.080000)" ### P(otential) E(vapo)T(ranspiration)