[GRASS-user] [GRASSLIST:1178] how to use cloud cover with r.sun?

GRASS GRASS user
1

Hello good people,

Can anyone suggest how I could use cloud cover (percentage) with r.sun, to generate global radiation maps which account for cloud cover?

I dont believe I can get the data to calibrate the a and b coefficients as mentioned in the manual. This is why I am hoping there is a simple equation that allows me to use cloud cover percentage directly…I have the cloud cover maps.

FYI, I am using Grass 6.0.1 on Ubuntu 5.1.

And I am attempting to generate monthly Potential Evapotranspiration maps from 1960 to 2002 for the Indian subcontinent.

cheers!
vishal

2

On Monday 21 August 2006 22:00, Vishal Mehta wrote:

Hello good people,

Can anyone suggest how I could use cloud cover (percentage) with r.sun, to
generate global radiation maps which account for cloud cover?

I dont believe I can get the data to calibrate the a and b coefficients as
mentioned in the manual. This is why I am hoping there is a simple equation
that allows me to use cloud cover percentage directly..I have the cloud
cover maps.

FYI, I am using Grass 6.0.1 on Ubuntu 5.1.

And I am attempting to generate monthly Potential Evapotranspiration maps
from 1960 to 2002 for the Indian subcontinent.

cheers!
vishal

Interesting questions. I have found that for cases in which cloud cover data
is missing, a best-estimate of the Linke Turbidity Factor can make large
improvements in the modeled daily summed radiation. see attached comparison
of radiation as measured by a WX station.

From the manual it looks like estimations of the _real_ atmospheric condition
must be described by two values: coefbh and coefdh :

Besides clear-sky radiations, user can compute a real-sky radiation (beam,
diffuse) using coefbh and coefdh input raster maps defining the fraction of
the respective clear-sky radiations reduced by atmospheric factors (e.g.
cloudiness). The value is between 0-1. Usually these coefficients can be
obtained from a long-terms meteorological measurements.

Good luck!

--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341

(attachments)

2005-initial_results.png

3

Hi,
You can use clear-sky index map and multiply it by your r.sun output.
I usually calculate clear-sky index as the ratio between the 'clear sky condition radiation (the r.sun)' and the 'observed radiation'
then I interpolate the point to retrive the map.

Actually I don't really know what is the influence of cloud cover on radiation (linear function between cloud cover and energy filtered???), probably you could calculate the reölation from observations, and then apply this with r.mapcalc.

Hope it helps,
Massimiliano

Dylan Beaudette wrote:

On Monday 21 August 2006 22:00, Vishal Mehta wrote:
  

Hello good people,

Can anyone suggest how I could use cloud cover (percentage) with r.sun, to
generate global radiation maps which account for cloud cover?

I dont believe I can get the data to calibrate the a and b coefficients as
mentioned in the manual. This is why I am hoping there is a simple equation
that allows me to use cloud cover percentage directly..I have the cloud
cover maps.

FYI, I am using Grass 6.0.1 on Ubuntu 5.1.

And I am attempting to generate monthly Potential Evapotranspiration maps
from 1960 to 2002 for the Indian subcontinent.

cheers!
vishal
    
Interesting questions. I have found that for cases in which cloud cover data is missing, a best-estimate of the Linke Turbidity Factor can make large improvements in the modeled daily summed radiation. see attached comparison of radiation as measured by a WX station.

From the manual it looks like estimations of the _real_ atmospheric condition must be described by two values: coefbh and coefdh :

Besides clear-sky radiations, user can compute a real-sky radiation (beam, diffuse) using coefbh and coefdh input raster maps defining the fraction of the respective clear-sky radiations reduced by atmospheric factors (e.g. cloudiness). The value is between 0-1. Usually these coefficients can be obtained from a long-terms meteorological measurements.

Good luck!

------------------------------------------------------------------------

------------------------------------------------------------------------

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser
  
--

Dr. Eng. Massimiliano Cannata
Scuola Universitaria Professionale della Svizzera Italiana
Istituto Scienze della Terra
Via Trevano, c.p. 72
CH-6952 Canobbio-Lugano
Tel: +41 (0)58 666 62 14
Fax +41 (0)58 666 62 09

4

Massimiliano Cannata wrote:

Actually I don't really know what is the influence of cloud cover on
radiation (linear function between cloud cover and energy
filtered???),

on a cloudy day the direct sunlight is reduced, and the diffuse in
increased. e.g. in the fjords I study there are many places which
never get direct sunlight. We have deployed subsurface PAR light-loggers
in paired sunny/shaded sites and the shaded side actually gets more
light on cloudy days.

Hamish

5

Hi all,
Right now i’ll go ahead with clear sky algorithms, and i’ll order observed data and try to come up with a real-sky correction based on cloud cover, r.sun output and observed data…
But i noticed a major problem as i ran r.sun with the shadowing effect (-s) flag:
Much of the centre of India ends up with 0 beam radiation!!
See attached images- which are beam radiation output for the same day, with and without shadow. yellow color is zero beam radiation.

THE REASON FOR THIS (I think)
I used the GTOPO30 (approx. 1km ORIGINAL res) elevation and derived slope and aspect as input to r.sun. HOWEVER, i never used the original GTOPO30 resolution- I just brought it into my current region, which is based on my base datasets of meterology from CRU at 0.5deg resoluton (approx. 50km? resolution), using r.in.bin directly into my coarse region (0.5 deg).
I think as a result, the shadowing effect in r.sun, which takes the max altitude from 4 nearest cells, is reaching across 50km(!) Hence with the Himalayas in the north and Western Ghats mountains in the west, most of central India gets 0 beam radiation!
At this point, I would like some help in how i go about this.

  1. How Can i import the GTOPO30 into my current mapset while at the same time creating a new region (say gtoporegion)?

  2. Then i think i should calculate my slopes and aspects at the orginal GTOP30s resolution (30arcseconds) in gtoporegion

  3. Then I should calculate all the r.sun maps in this gtoporegion

  4. Finally I guess I can take these fine resolution solar maps into my 0.5degree resolution region using r.resample.rst, and continue …??

Any suggestions appreciated! I’m not very familiar with having different resolutions/regions in the same mapset (have read the documentation).

thanks,
vishal

On 8/23/06, Hamish <hamish_nospam@yahoo.com> wrote:

Massimiliano Cannata wrote:

Actually I don’t really know what is the influence of cloud cover on
radiation (linear function between cloud cover and energy
filtered???),

on a cloudy day the direct sunlight is reduced, and the diffuse in
increased. e.g. in the fjords I study there are many places which
never get direct sunlight. We have deployed subsurface PAR light-loggers
in paired sunny/shaded sites and the shaded side actually gets more
light on cloudy days.

Hamish

grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

(attachments)

withshadow.png
withoutshadow.png

6

Vishal Mehta wrote:

But i noticed a major problem as i ran r.sun with the shadowing effect
(-s) flag:
Much of the centre of India ends up with 0 beam radiation!!
See attached images- which are beam radiation output for the same day,
with and without shadow. yellow color is zero beam radiation.

THE REASON FOR THIS (I think)
I used the GTOPO30 (approx. 1km ORIGINAL res) elevation and derived
slope and aspect as input to r.sun. HOWEVER, i never used the original
GTOPO30 resolution- I just brought it into my current region, which is
based on my base datasets of meterology from CRU at 0.5deg resoluton
(approx. 50km? resolution), using r.in.bin directly into my coarse
region (0.5 deg). I think as a result, the shadowing effect in r.sun,
which takes the max altitude from 4 nearest cells, is reaching across
50km(!) Hence with the Himalayas in the north and Western Ghats
mountains in the west, most of central India gets 0 beam radiation!

I think you are correct, elevation in meters is probably being treaded
as degrees, so you get an unwanted 1852*60 vertical exaggeration.
Sounds plausible at least.

Solutions:

a) reproject your DEM into a meter-based projection (I couldn't tell you
which one), then reproject the results back to lat/lon.

b) rescale your DEM into "degrees". e.g.:

r.mapcalc elevmap_scaled="elevmap / 1852*60"
  [maybe figure in cos(mean_lat)]

That will probably mess up the linke turbidity functions which I believe
takes physical elevation into account. (but I could be wrong)

Hamish

7

Hamish,

You may be right that elevation in meters is being used as degrees (my region is lat long, elevation in meters)…because since my last email I have run the r.sun module using a new region based on the finer GTOPO30 resolution. With the shadow option on, i am getting even more zeros- almost all zeros!!

But is it really possible that the r.sun module (i’m using Grass 6.0.1 on Ubuntu Breezy) would treat elevation in degrees??

And can you please tell me where the 1852*60 comes from?
I guess i would have to creat a new location in UTM coordinates and then do the projecting back and forth you suggest…

Jose, I am not using a particular model for PET, rather I am using a bunch of equations mostly from FAO-56 (just google it, the book is online); met. data from the Climate Research Unit, UK.
And there’s no easy way to do it, at least based on my research. To give just 2 examples,

  1. r.sun module gives solar radiation(clear sky), and you need observed spatio-temporal data (difficult to get) on either clear sky index or the coeffd, coefb parameters to get real sky radiation;
  2. The FAO 56 manual does have some simple options. but these involve working with extra terrestrial radiation- i dont think r.sun gives that (its internally calculated i believe).

anyway…it might take me 2-3 months to produce some final output…

cheers,
vishal

On 8/24/06, Hamish <hamish_nospam@yahoo.com> wrote:

I think you are correct, elevation in meters is probably being treaded
as degrees, so you get an unwanted 1852*60 vertical exaggeration.
Sounds plausible at least.

Solutions:

a) reproject your DEM into a meter-based projection (I couldn’t tell you
which one), then reproject the results back to lat/lon.

b) rescale your DEM into “degrees”. e.g.:

r.mapcalc elevmap_scaled=“elevmap / 1852*60”
[maybe figure in cos(mean_lat)]

That will probably mess up the linke turbidity functions which I believe
takes physical elevation into account. (but I could be wrong)

Hamish

8

On Thu, Aug 24, 2006 at 04:40:33PM +0530, Vishal Mehta wrote:

   Hamish,

   You may be right that elevation in meters is being used as degrees (my
   region is lat long, elevation in meters)...

As far as I know the latlong usage is disabled in r.topidx
exactly for this reason: r.topidx cannot deal with latlong
locations and elevation being in meters.

Maybe the same applies to r.sun? I have added Jaro in cc,
probably he can send us a comment.

Markus

PS: Jaro: the thread is here:
    http://grass.itc.it/pipermail/grassuser/2006-August/thread.html#35863

9

Vishal Mehta wrote:

You may be right that elevation in meters is being used as degrees (my
region is lat long, elevation in meters)...because since my last email I
have run the r.sun module using a new region based on the finer GTOPO30
resolution. With the shadow option on, i am getting even more zeros- almost
all zeros!!

But is it really possible that the r.sun module (i'm using Grass 6.0.1 on
Ubuntu Breezy) would treat elevation in degrees??

There is a tendency to assume that all distances are in the same
(unspecified) units, so if two elevations differ by one unit and are
one unit apart on the ground, the slope is 45 degrees.

If both coordinates and elevations are in metres, there's no problem.
If they're both in feet, again no problem. If one is in feet and the
other in metres, the results will be wrong. If the coordinates are in
degrees, the results are likely to be very wrong.

This problem can be solved to an extent by using the value from the
PROJ_UNITS file to convert ground distances to metres (although that
requires that the PROJ_UNITS file actually contains the correct
value).

--
Glynn Clements <glynn@gclements.plus.com>

10

On Thursday 24 August 2006 09:29, Glynn Clements wrote:

Vishal Mehta wrote:
> You may be right that elevation in meters is being used as degrees (my
> region is lat long, elevation in meters)...because since my last email I
> have run the r.sun module using a new region based on the finer GTOPO30
> resolution. With the shadow option on, i am getting even more zeros-
> almost all zeros!!
>
> But is it really possible that the r.sun module (i'm using Grass 6.0.1 on
> Ubuntu Breezy) would treat elevation in degrees??

There is a tendency to assume that all distances are in the same
(unspecified) units, so if two elevations differ by one unit and are
one unit apart on the ground, the slope is 45 degrees.

If both coordinates and elevations are in metres, there's no problem.
If they're both in feet, again no problem. If one is in feet and the
other in metres, the results will be wrong. If the coordinates are in
degrees, the results are likely to be very wrong.

This problem can be solved to an extent by using the value from the
PROJ_UNITS file to convert ground distances to metres (although that
requires that the PROJ_UNITS file actually contains the correct
value).

Perhaps I have missed something here, but why not project you latlong data to
some coordinate system with x,y,z units in meters ?

--
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341

11

On Thu, Aug 24, 2006 at 10:02:23AM -0700, Dylan Beaudette wrote:
...

Perhaps I have missed something here, but why not project you latlong data to
some coordinate system with x,y,z units in meters ?

This should be a projection covering the entire subcontinent
in a reasonable way - what are the official projection(s) for
that area?

Markus

12

But is it really possible that the r.sun module (i'm using Grass 6.0.1
on Ubuntu Breezy) would treat elevation in degrees??

Yes, it is possible. Maybe even likely.

And can you please tell me where the 1852*60 comes from?

There are 1852 meters in a nautical mile. There are 60 nautical miles in
a degree of latitude. elev_in_meters / (1850.0 * 60) = elevation in deg

The width of a deg longitude is equal to width_lat*cos(lat).

If rescaling the elevation into degress be sure to /(1852.0*60). This
ensures that the result will be a floating point value (as opposed to
/(1852*60) ).

I guess i would have to creat a new location in UTM coordinates and
then do the projecting back and forth you suggest..

I think this is the saner solution, but I expect that a single UTM zone
is not an appropriate projection for working with an area the size of
India.

Markus wrote:

This should be a projection covering the entire subcontinent
in a reasonable way - what are the official projection(s) for
that area?

If you have problems finding one, you might ask for help on the PROJ.4
mailing list.

Hamish

13

Thanks everyone,
a few thoughts-

  1. Neither the r.sun manual, nor the r.slope.aspect manual, states that the x,y and elevation units should be the same. In fact both manuals say the elevation should be in meters! r.sun talks about how if the maps are projected then it does’nt matter- but it does not talk of use with lat-long system(nor does r.slope.aspect).

do these manuals then just assume that we are using a projection where the x,y are in meters??

  1. Hamish and Markus are right - my study area covers several utm zones.

  2. At this point i’ll try the dem, slope, aspect and r.sun in Lambert conformal conic projection with everest datum;
    or if that does not work, i’ll try Hamish’ idea of rescaling the elevation to degrees.

I still cant get over the idea that GRASS cant handle correctly slope generation in lat-long coordinates with elevation in meters!

thanks again,
vishal

On 8/25/06, Hamish <hamish_nospam@yahoo.com > wrote:

But is it really possible that the r.sun module (i’m using Grass 6.0.1
on Ubuntu Breezy) would treat elevation in degrees??

Yes, it is possible. Maybe even likely.

And can you please tell me where the 1852*60 comes from?

There are 1852 meters in a nautical mile. There are 60 nautical miles in
a degree of latitude. elev_in_meters / (1850.0 * 60) = elevation in deg

The width of a deg longitude is equal to width_lat*cos(lat).

If rescaling the elevation into degress be sure to /(1852.060). This
ensures that the result will be a floating point value (as opposed to
/(185260) ).

I guess i would have to creat a new location in UTM coordinates and
then do the projecting back and forth you suggest…

I think this is the saner solution, but I expect that a single UTM zone
is not an appropriate projection for working with an area the size of
India.

Markus wrote:

This should be a projection covering the entire subcontinent
in a reasonable way - what are the official projection(s) for
that area?

If you have problems finding one, you might ask for help on the PROJ.4
mailing list.

Hamish

14

On 8/24/06, Glynn Clements <glynn@gclements.plus.com> wrote:

This problem can be solved to an extent by using the value from the
PROJ_UNITS file to convert ground distances to metres (although that
requires that the PROJ_UNITS file actually contains the correct
value).

I did’nt fully understand this comment Glynn. Could you please elaborate this idea?
FYI,
GRASS 6.0.1 (cru3):~/CRU > g.proj -w
GEOGCS[“wgs84”,
DATUM[“WGS_1984”,
SPHEROID[“wgs84”,6378137,298.257223563],
TOWGS84[0.000,0.000,0.000]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433]]

And my Proj_units file in the PERMANENT mapset is:
unit: degree
units: degrees
meters: 1.0

I’m not sure what the last line of g.proj -w means, what does unit there refer to,
as opposed to the proj_units ‘meters’ field?

thanks,
vishal

15

Vishal Mehta wrote:

> This problem can be solved to an extent by using the value from the
> PROJ_UNITS file to convert ground distances to metres (although that
> requires that the PROJ_UNITS file actually contains the correct
> value).

I did'nt fully understand this comment Glynn. Could you please elaborate
this idea?
FYI,
GRASS 6.0.1 (cru3):~/CRU > g.proj -w
GEOGCS["wgs84",
    DATUM["WGS_1984",
        SPHEROID["wgs84",6378137,298.257223563],
        TOWGS84[0.000,0.000,0.000]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

And my Proj_units file in the PERMANENT mapset is:

unit: degree
units: degrees
meters: 1.0

The "meters" value in the PROJ_UNITS file is supposed to be a scale
factor which converts map units to metres (i.e. 1.0 for metres, 0.3048
for feet). For degrees, it should be roughly 111300 * cos(latitude).

r.slope.aspect uses that value to convert map units to metres; if it's
set correctly (and the zfactor= option is set correctly if your
elevations aren't in metres), your slope map should be correct. There
might still be problems with r.sun itself, though.

I'm not sure what the last line of g.proj -w means, what does unit there
refer to, as opposed to the proj_units 'meters' field?

0.0174532925199433 is pi/180, i.e. it represents the size of the
(angular) units in radians.

--
Glynn Clements <glynn@gclements.plus.com>

16

Glynn Clements wrote:

> > This problem can be solved to an extent by using the value from the
> > PROJ_UNITS file to convert ground distances to metres (although that
> > requires that the PROJ_UNITS file actually contains the correct
> > value).
>
> I did'nt fully understand this comment Glynn. Could you please elaborate
> this idea?
> FYI,
> GRASS 6.0.1 (cru3):~/CRU > g.proj -w
> GEOGCS["wgs84",
> DATUM["WGS_1984",
> SPHEROID["wgs84",6378137,298.257223563],
> TOWGS84[0.000,0.000,0.000]],
> PRIMEM["Greenwich",0],
> UNIT["degree",0.0174532925199433]]
>
> And my Proj_units file in the PERMANENT mapset is:

> unit: degree
> units: degrees
> meters: 1.0

The "meters" value in the PROJ_UNITS file is supposed to be a scale
factor which converts map units to metres (i.e. 1.0 for metres, 0.3048
for feet). For degrees, it should be roughly 111300 * cos(latitude).

r.slope.aspect uses that value to convert map units to metres; if it's
set correctly (and the zfactor= option is set correctly if your
elevations aren't in metres), your slope map should be correct. There
might still be problems with r.sun itself, though.

Ignore that last part; r.slope.aspect actually ignores that value, and
instead uses G_distance(), which has a special case for lat/lon.

--
Glynn Clements <glynn@gclements.plus.com>

17

so should i change the Proj_units file to:

unit: degree
units: degrees
meters: 111300 * cos(y())

-vishal

On 8/25/06, Glynn Clements <glynn@gclements.plus.com> wrote:

Vishal Mehta wrote:

This problem can be solved to an extent by using the value from the
PROJ_UNITS file to convert ground distances to metres (although that
requires that the PROJ_UNITS file actually contains the correct
value).

I did’nt fully understand this comment Glynn. Could you please elaborate
this idea?
FYI,
GRASS 6.0.1 (cru3):~/CRU > g.proj -w
GEOGCS[“wgs84”,
DATUM[“WGS_1984”,
SPHEROID[“wgs84”,6378137,298.257223563],
TOWGS84[0.000,0.000,0.000]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”, 0.0174532925199433]]

And my Proj_units file in the PERMANENT mapset is:

unit: degree
units: degrees
meters: 1.0

The “meters” value in the PROJ_UNITS file is supposed to be a scale
factor which converts map units to metres (i.e. 1.0 for metres, 0.3048
for feet). For degrees, it should be roughly 111300 * cos(latitude).

r.slope.aspect uses that value to convert map units to metres; if it’s
set correctly (and the zfactor= option is set correctly if your
elevations aren’t in metres), your slope map should be correct. There
might still be problems with r.sun itself, though.

I’m not sure what the last line of g.proj -w means, what does unit there
refer to, as opposed to the proj_units ‘meters’ field?

0.0174532925199433 is pi/180, i.e. it represents the size of the
(angular) units in radians.

–
Glynn Clements < glynn@gclements.plus.com>

18

Vishal Mehta wrote:

so should i change the Proj_units file to:

unit: degree
units: degrees
meters: 111300 * cos(y())

It won't matter in this case; both r.slope.aspect and r.sun handle
lat/lon locations separately.

--
Glynn Clements <glynn@gclements.plus.com>

19

Hi all,

the idea of converting elevation to degrees did’nt work out.

i managed to project my latlong dem (30min, approx 960m resolution) to lambert conformal conic, so elevations, and coordinates are both in meters.

then i ran r.sun with a fixed albedo of 0.23

r.sun elevin=gtopo30 aspin=gtopoasp slopein=gtoposlope alb=0.23 beam_rad=beam1 diff_rad=diffuse1 refl_rad=reflect1 day=15

i got some output that looks right for beam and diffuse radiation: however, the reflected radiation map is almost all zero. there are some small positive values near mountainous regions (~50 Whr/(m2-day)).

since the albedo was 0.23, and the average global radiation i get is about 5000 Whr/(m2-day), i cant understand why the reflected radiation (Ri)is all zero. The manual says something about a factor related to ‘fraction of ground viewed by an inclined surface’ rg(gN).

i.e. Ri = albedo* global radiation * rg(gN)

so I assume then that this rg(gN) factor must be causing the problem…

does anyone have any idea why this may be happening and what to do about it?

also, i wanted to report that what i’m getting for beam, diffuse and reflected output now with lambert projections and the original 30min (960m) resolution float dem,
is almost the SAME as in the beginning of this saga, when i was using latlong at coarse 0.5 degree resolution and integer dem! The difference may only be in the output when using the -s flag for accounting for shading.
However, the -s flag option took some half hour to run at the 960m resolution…

There are many other things that happened in between that were funky, but i guess i should save you all from that stuff:)

cheers,
vishal

On 8/28/06, Glynn Clements <glynn@gclements.plus.com> wrote:

Vishal Mehta wrote:

so should i change the Proj_units file to:

unit: degree
units: degrees
meters: 111300 * cos(y())

It won’t matter in this case; both r.slope.aspect and r.sun handle
lat/lon locations separately.

–
Glynn Clements <glynn@gclements.plus.com>