[GRASS-SVN] r38772 - grass/trunk/raster/r.sun
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Aug 17 09:20:59 EDT 2009
Author: hamish
Date: 2009-08-17 09:20:58 -0400 (Mon, 17 Aug 2009)
New Revision: 38772
Modified:
grass/trunk/raster/r.sun/main.c
grass/trunk/raster/r.sun/r.sun.html
grass/trunk/raster/r.sun/rsunlib.c
Log:
trac #498: (merge from devbr6)
lat= is no longer usable- chop it out;
make latin= and longin= work, but they give us only ~1.3% speed boost);
give some variables more understandable names
Modified: grass/trunk/raster/r.sun/main.c
===================================================================
--- grass/trunk/raster/r.sun/main.c 2009-08-17 13:19:15 UTC (rev 38771)
+++ grass/trunk/raster/r.sun/main.c 2009-08-17 13:20:58 UTC (rev 38772)
@@ -126,11 +126,11 @@
int varCount_global = 0;
int bitCount_global = 0;
int arrayNumInt = 1;
-float **z = NULL, **o = NULL, **s = NULL, **li = NULL, **a = NULL, **la =
- NULL, **longitArray, **cbhr = NULL, **cdhr = NULL;
+float **z = NULL, **o = NULL, **s = NULL, **li = NULL, **a = NULL,
+ **latitudeArray = NULL, **longitudeArray, **cbhr = NULL, **cdhr = NULL;
double op, dp;
double invstepx, invstepy;
-double sr_min = 24., sr_max = 0., ss_min = 24., ss_max = 0.;
+double sunrise_min = 24., sunrise_max = 0., sunset_min = 24., sunset_max = 0.;
float **lumcl, **beam, **insol, **diff, **refl, **globrad;
unsigned char *horizonarray = NULL;
@@ -141,10 +141,10 @@
*/
double xmin, xmax, ymin, ymax;
double declin, step, dist;
-double li_max = 0., li_min = 100., al_max = 0., al_min = 1.0, la_max = -90.,
- la_min = 90.;
+double linke_max = 0., linke_min = 100., albedo_max = 0., albedo_min = 1.0,
+ lat_max = -90., lat_min = 90.;
double offsetx = 0.5, offsety = 0.5;
-char *tt, *lt;
+char *ttime;
/*
* double slope;
@@ -208,7 +208,7 @@
struct
{
struct Option *elevin, *aspin, *aspect, *slopein, *slope, *linkein,
- *lin, *albedo, *longin, *alb, *latin, *lat, *coefbh, *coefdh,
+ *lin, *albedo, *longin, *alb, *latin, *coefbh, *coefdh,
*incidout, *beam_rad, *insol_time, *diff_rad, *refl_rad,
*glob_rad, *day, *step, *declin, *ltime, *dist, *horizon,
*horizonstep, *numPartitions, *civilTime;
@@ -328,26 +328,16 @@
parm.latin->required = NO;
parm.latin->gisprompt = "old,cell,raster";
parm.latin->description =
- _("Name of the latitudes input raster map [decimal degrees]");
+ _("Name of input raster map containing latitudes [decimal degrees]");
parm.latin->guisection = _("Input_options");
- if (parm.latin->answer == NULL) {
- parm.lat = G_define_option();
- parm.lat->key = "lat";
- parm.lat->type = TYPE_DOUBLE;
- parm.lat->required = NO;
- parm.lat->description =
- _("A single value of latitude [decimal degrees]");
- parm.lat->guisection = _("Input_options");
- }
-
parm.longin = G_define_option();
parm.longin->key = "longin";
parm.longin->type = TYPE_STRING;
parm.longin->required = NO;
parm.longin->gisprompt = "old,cell,raster";
parm.longin->description =
- _("Name of the longitude input raster map [decimal degrees]");
+ _("Name of input raster map containing longitudes [decimal degrees]");
parm.longin->guisection = _("Input_options");
parm.coefbh = G_define_option();
@@ -553,20 +543,20 @@
linkein = parm.linkein->answer;
albedo = parm.albedo->answer;
latin = parm.latin->answer;
+ longin = parm.longin->answer;
civiltime = parm.civilTime->answer;
if (civiltime != NULL) {
setUseCivilTime(TRUE);
- longin = parm.longin->answer;
- if (longin == NULL) {
- G_fatal_error(_("You must give the longitude raster if you use civil time"));
+ if (longin == NULL)
+ G_fatal_error(
+ _("You must give the longitude raster if you use civil time"));
- }
-
if (sscanf(parm.civilTime->answer, "%lf", &civilTime) != 1)
G_fatal_error(_("Error reading civil time zone value"));
+
if (civilTime < -24. || civilTime > 24.)
G_fatal_error(_("Invalid civil time zone value"));
@@ -614,7 +604,7 @@
G_fatal_error(_("If you use the horizon option you must also set the 'horizonstep' parameter."));
}
- tt = parm.ltime->answer;
+ ttime = parm.ltime->answer;
if (parm.ltime->answer != NULL) {
if (insol_time != NULL)
G_fatal_error(_("Time and insol_time are incompatible options"));
@@ -651,15 +641,6 @@
sscanf(parm.lin->answer, "%lf", &singleLinke);
if (parm.albedo->answer == NULL)
sscanf(parm.alb->answer, "%lf", &singleAlbedo);
- lt = parm.lat->answer;
- /*
- * if (parm.lat->answer != NULL)
- * sscanf(parm.lat->answer, "%lf", &latitude);
- */
- /* HB 6/2008: why is the above commented out? instead of sscanf, maybe nicer to use:
- G_scan_lat(parm.lat->answer, &latitude);
- MN 2/2009: latitude doesn't exist! also G_scan_lat() does not exist in GRASS
- */
if (parm.slopein->answer == NULL)
sscanf(parm.slope->answer, "%lf", &singleSlope);
@@ -693,7 +674,9 @@
* shadows pre-calculated, there is no problem. */
if (saveMemory && useShadow() && (!useHorizonData()))
- G_fatal_error(_("If you want to save memory and to use shadows, you must use pre-calculated horizons."));
+ G_fatal_error(
+ _("If you want to save memory and to use shadows, "
+ "you must use pre-calculated horizons."));
if (parm.declin->answer == NULL)
declination = com_declin(day);
@@ -702,12 +685,7 @@
declination = -declin;
}
-
- /*
- * if (lt != NULL)
- * latitude = -latitude * deg2rad;
- */
- if (tt != 0) {
+ if (ttime != 0) {
/* Shadow for just one time during the day */
if (horizon == NULL) {
arrayNumInt = 1;
@@ -723,7 +701,7 @@
}
}
- if (tt != NULL) {
+ if (ttime != NULL) {
tim = (timo - 12) * 15;
/* converting to degrees */
@@ -735,19 +713,22 @@
}
/* Set up parameters for projection to lat/long if necessary */
-
- if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {
+ /* we need to do this even if latin= and longin= were given because
+ rsunlib's com_par() needs iproj and oproj */
+ if ((G_projection() != PROJECTION_LL)) {
struct Key_Value *in_proj_info, *in_unit_info;
if ((in_proj_info = G_get_projinfo()) == NULL)
- G_fatal_error
- (_("Can't get projection info of current location: please set latitude via 'lat' or 'latin' option!"));
+ G_fatal_error(
+ _("Can't get projection info of current location"));
if ((in_unit_info = G_get_projunits()) == NULL)
- G_fatal_error(_("Can't get projection units of current location"));
+ G_fatal_error(
+ _("Can't get projection units of current location"));
if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
- G_fatal_error(_("Can't get projection key values of current location"));
+ G_fatal_error(
+ _("Can't get projection key values of current location"));
G_free_key_value(in_proj_info);
G_free_key_value(in_unit_info);
@@ -760,8 +741,15 @@
G_fatal_error(_("Unable to set up lat/long projection parameters"));
}
+ if ((latin != NULL || longin != NULL) && (G_projection() == PROJECTION_LL))
+ G_warning(_("latin and longin raster maps have no effect when in a Lat/Lon location"));
+ /* true about longin= when civiltime is used? */
+ if ((latin == NULL && longin != NULL) || (latin != NULL && longin == NULL))
+ G_fatal_error(_("Both latin and longin raster maps must be given, or neither"));
+
/**********end of parser - ******************************/
+
if ((G_projection() == PROJECTION_LL))
ll_correction = TRUE;
@@ -785,8 +773,8 @@
FCELL *rast1 = NULL, *rast2 = NULL;
static FCELL **horizonbuf;
unsigned char *horizonpointer;
- int fd1 = -1, fd2 = -1, fd3 = -1, fd4 = -1, fd5 = -1, fd6 = -1, fd7 =
- -1, row, row_rev;
+ int fd1 = -1, fd2 = -1, fd3 = -1, fd4 = -1, fd5 = -1, fd6 = -1,
+ fd7 = -1, row, row_rev;
static int *fd_shad;
int fr1 = -1, fr2 = -1;
int l, i, j;
@@ -877,10 +865,10 @@
if (latin != NULL) {
cell6 = Rast_allocate_f_buf();
- if (la == NULL) {
- la = (float **)G_malloc(sizeof(float *) * (numRows));
+ if (latitudeArray == NULL) {
+ latitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
for (l = 0; l < numRows; l++)
- la[l] = (float *)G_malloc(sizeof(float) * (n));
+ latitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
}
if ((mapset = G_find_raster2(latin, "")) == NULL)
G_fatal_error(_("Raster map <%s> not found"), latin);
@@ -889,9 +877,9 @@
if (longin != NULL) {
cell7 = Rast_allocate_f_buf();
- longitArray = (float **)G_malloc(sizeof(float *) * (numRows));
+ longitudeArray = (float **)G_malloc(sizeof(float *) * (numRows));
for (l = 0; l < numRows; l++)
- longitArray[l] = (float *)G_malloc(sizeof(float) * (n));
+ longitudeArray[l] = (float *)G_malloc(sizeof(float) * (n));
if ((mapset = G_find_raster2(longin, "")) == NULL)
G_fatal_error(_("Raster map <%s> not found"), longin);
@@ -932,7 +920,7 @@
fd_shad = (int *)G_malloc(sizeof(int) * arrayNumInt);
}
/*
- * if(tt != NULL)
+ * if(ttime != NULL)
* {
*
* horizonbuf[0]=Rast_allocate_f_buf();
@@ -1041,16 +1029,16 @@
if (latin != NULL) {
if (!Rast_is_f_null_value(cell6 + j))
- la[rowrevoffset][j] = (float)cell6[j];
+ latitudeArray[rowrevoffset][j] = (float)cell6[j];
else
- la[rowrevoffset][j] = UNDEFZ;
+ latitudeArray[rowrevoffset][j] = UNDEFZ;
}
if (longin != NULL) {
if (!Rast_is_f_null_value(cell7 + j))
- longitArray[rowrevoffset][j] = (float)cell7[j];
+ longitudeArray[rowrevoffset][j] = (float)cell7[j];
else
- longitArray[rowrevoffset][j] = UNDEFZ;
+ longitudeArray[rowrevoffset][j] = UNDEFZ;
}
if (coefbh != NULL) {
@@ -1141,7 +1129,7 @@
z[i][j] = UNDEFZ;
if (albedo != NULL && a[i][j] == UNDEFZ)
z[i][j] = UNDEFZ;
- if (latin != NULL && la[i][j] == UNDEFZ)
+ if (latin != NULL && latitudeArray[i][j] == UNDEFZ)
z[i][j] = UNDEFZ;
if (coefbh != NULL && cbhr[i][j] == UNDEFZ)
z[i][j] = UNDEFZ;
@@ -1335,7 +1323,7 @@
com_par(sunGeom, sunVarGeom, gridGeom, latitude, longitude);
- if (tt != NULL) { /*irradiance */
+ if (ttime != NULL) { /*irradiance */
s0 = lumcline2(sunGeom, sunVarGeom, sunSlopeGeom, gridGeom,
horizonpointer);
@@ -1810,7 +1798,8 @@
for (i = 0; i < n; i++) {
if (useCivilTime()) {
- longitTime = -longitArray[arrayOffset][i] / 15.;
+ /* sun travels 15deg per hour, so 1 TZ every 15 deg and 15 TZs * 24hrs = 360deg */
+ longitTime = -longitudeArray[arrayOffset][i] / 15.;
}
gridGeom.xg0 = gridGeom.xx0 = (double)i *gridGeom.stepx;
@@ -1841,42 +1830,49 @@
}
if (linkein != NULL) {
sunRadVar.linke = li[arrayOffset][i];
- li_max = AMAX1(li_max, sunRadVar.linke);
- li_min = AMIN1(li_min, sunRadVar.linke);
+ linke_max = AMAX1(linke_max, sunRadVar.linke);
+ linke_min = AMIN1(linke_min, sunRadVar.linke);
}
if (albedo != NULL) {
sunRadVar.alb = a[arrayOffset][i];
- al_max = AMAX1(al_max, sunRadVar.alb);
- al_min = AMIN1(al_min, sunRadVar.alb);
+ albedo_max = AMAX1(albedo_max, sunRadVar.alb);
+ albedo_min = AMIN1(albedo_min, sunRadVar.alb);
}
if (latin != NULL) {
- latitude = la[arrayOffset][i];
- la_max = AMAX1(la_max, latitude);
- la_min = AMIN1(la_min, latitude);
+ latitude = latitudeArray[arrayOffset][i];
+ lat_max = AMAX1(lat_max, latitude);
+ lat_min = AMIN1(lat_min, latitude);
latitude *= deg2rad;
}
- /* MN 2/2009: should it be??
- if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {
- */
+ if (longin != NULL) {
+ longitude = longitudeArray[arrayOffset][i];
+ /* lon_max = AMAX1(lon_max, longitude); */
+ /* lon_min = AMIN1(lon_min, longitude); */
+ longitude *= deg2rad;
+ }
+
if ((G_projection() != PROJECTION_LL)) {
- longitude = gridGeom.xp;
- latitude = gridGeom.yp;
+ if (latin == NULL || longin == NULL) {
+ /* if either is missing we have to calc both from current projection */
+ longitude = gridGeom.xp;
+ latitude = gridGeom.yp;
- if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
- G_fatal_error("Error in pj_do_proj");
- }
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0) {
+ G_fatal_error("Error in pj_do_proj");
+ }
- la_max = AMAX1(la_max, latitude);
- la_min = AMIN1(la_min, latitude);
- latitude *= deg2rad;
- longitude *= deg2rad;
+ lat_max = AMAX1(lat_max, latitude);
+ lat_min = AMIN1(lat_min, latitude);
+ latitude *= deg2rad;
+ longitude *= deg2rad;
+ }
}
else { /* ll projection */
latitude = gridGeom.yp;
longitude = gridGeom.xp;
- la_max = AMAX1(la_max, latitude);
- la_min = AMIN1(la_min, latitude);
+ lat_max = AMAX1(lat_max, latitude);
+ lat_min = AMIN1(lat_min, latitude);
latitude *= deg2rad;
longitude *= deg2rad;
}
@@ -1892,7 +1888,7 @@
cos_v = cos(M_PI / 2 + sunSlopeGeom.aspect);
sin_v = sin(M_PI / 2 + sunSlopeGeom.aspect);
- if (tt != NULL)
+ if (ttime != NULL)
sunGeom.timeAngle = tim;
gridGeom.sinlat = sin(-latitude);
@@ -1912,10 +1908,10 @@
if ((incidout != NULL) || someRadiation) {
com_par_const(longitTime, &sunGeom, &gridGeom);
- sr_min = AMIN1(sr_min, sunGeom.sunrise_time);
- sr_max = AMAX1(sr_max, sunGeom.sunrise_time);
- ss_min = AMIN1(ss_min, sunGeom.sunset_time);
- ss_max = AMAX1(ss_max, sunGeom.sunset_time);
+ sunrise_min = AMIN1(sunrise_min, sunGeom.sunrise_time);
+ sunrise_max = AMAX1(sunrise_max, sunGeom.sunrise_time);
+ sunset_min = AMIN1(sunset_min, sunGeom.sunset_time);
+ sunset_max = AMAX1(sunset_max, sunGeom.sunset_time);
}
@@ -1985,7 +1981,7 @@
day);
hist.edlinecnt = 2;
- if (tt != NULL) {
+ if (ttime != NULL) {
sprintf(hist.edhist[hist.edlinecnt],
" Local (solar) time (decimal hr.): %.4f", timo);
hist.edlinecnt++;
@@ -2000,17 +1996,12 @@
" Declination (rad): %f", -declination);
hist.edlinecnt += 3;
- if (lt != NULL)
- sprintf(hist.edhist[hist.edlinecnt],
- " Latitude (deg): %.4f",
- -latitude * rad2deg);
- else
- sprintf(hist.edhist[hist.edlinecnt],
- " Latitude min-max(deg): %.4f - %.4f",
- la_min, la_max);
+ sprintf(hist.edhist[hist.edlinecnt],
+ " Latitude min-max(deg): %.4f - %.4f",
+ lat_min, lat_max);
hist.edlinecnt++;
- if (tt != NULL) {
+ if (ttime != NULL) {
sprintf(hist.edhist[hist.edlinecnt],
" Sunrise time (hr.): %.2f",
sunGeom.sunrise_time);
@@ -2024,16 +2015,16 @@
else {
sprintf(hist.edhist[hist.edlinecnt],
" Sunrise time min-max (hr.): %.2f - %.2f",
- sr_min, sr_max);
+ sunrise_min, sunrise_max);
sprintf(hist.edhist[hist.edlinecnt + 1],
" Sunset time min-max (hr.): %.2f - %.2f",
- ss_min, ss_max);
+ sunset_min, sunset_max);
sprintf(hist.edhist[hist.edlinecnt + 2],
" Time step (hr.): %.4f", step);
}
hist.edlinecnt += 3;
- if (incidout != NULL || tt != NULL) {
+ if (incidout != NULL || ttime != NULL) {
sprintf(hist.edhist[hist.edlinecnt],
" Solar altitude (deg): %.4f",
sunVarGeom.solarAltitude * rad2deg);
@@ -2050,7 +2041,7 @@
else
sprintf(hist.edhist[hist.edlinecnt],
" Linke turbidity factor min-max: %.1f-%.1f",
- li_min, li_max);
+ linke_min, linke_max);
hist.edlinecnt++;
if (albedo == NULL)
@@ -2060,7 +2051,7 @@
else
sprintf(hist.edhist[hist.edlinecnt],
" Ground albedo min-max: %.3f-%.3f",
- al_min, al_max);
+ albedo_min, albedo_max);
hist.edlinecnt++;
sprintf(hist.edhist[hist.edlinecnt],
Modified: grass/trunk/raster/r.sun/r.sun.html
===================================================================
--- grass/trunk/raster/r.sun/r.sun.html 2009-08-17 13:19:15 UTC (rev 38771)
+++ grass/trunk/raster/r.sun/r.sun.html 2009-08-17 13:20:58 UTC (rev 38772)
@@ -44,7 +44,8 @@
parameters are saved in the resultant maps' history files, which may be viewed
with the <a href="r.info.html">r.info</a> command.
</p>
-<p>The solar incidence angle raster map <i>incidout</i> is computed specifying
+<p>
+The solar incidence angle raster map <i>incidout</i> is computed specifying
elevation raster map <i>elevin</i>, aspect raster map <i>aspin</i>, slope
steepness raster map <i>slopin,</i> given the day <i>day</i> and local time
<i>time</i>. There is no need to define latitude for locations with known
@@ -52,8 +53,7 @@
g.proj</a>
command). If you have undefined projection, (x,y) system, etc. then the latitude
can be defined explicitely for large areas by input raster map <i>latin</i>
- with interpolated latitude values or, for smaller areas, a single region
-latitude value <i>lat</i> can be used instead. All input raster maps must
+ with interpolated latitude values. All input raster maps must
be floating point (FCELL) raster maps. Null data in maps are excluded from
the computation (and also speeding-up the computation), so each output raster
map will contain null data in cells according to all input raster maps. The
@@ -63,7 +63,8 @@
where January 1 is day no.1 and December 31 is 365. Time <i>time</i> must
be a local (solar) time (i.e. NOT a zone time, e.g. GMT, CET) in decimal system,
e.g. 7.5 (= 7h 30m A.M.), 16.1 = 4h 6m P.M.. </p>
-<p>Setting the solar declination <i>declin</i> by user is an option to override
+<p>
+Setting the solar declination <i>declin</i> by user is an option to override
the value computed by the internal routine for the day of the year. The value
of geographical latitude can be set as a constant for the whole computed
region or, as an option, a grid representing spatially distributed values
@@ -71,15 +72,17 @@
with positive values for northern hemisphere and negative for southern one.
In similar principle the Linke turbidity factor (<i>linkein</i>, <i>lin</i>
) and ground albedo (<i>albedo</i>, <i>alb</i>) can be set. </p>
-<p>Besides clear-sky radiations, the user can compute a real-sky radiation (beam,
+<p>
+Besides clear-sky radiations, the user can compute a real-sky radiation (beam,
diffuse) using <i>coefbh</i> and <i>coefdh </i>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.</p>
-<p>The solar irradiation or irradiance raster maps <i>beam_rad</i>, <i>diff_rad</i>
-, <i>refl_rad</i> are computed for a given day <i>day,</i> latitude <i>lat
-(latin), </i>elevation <i>elevin</i>, slope <i>slopein</i> and aspect <i>
-aspin</i> raster maps. For convenience, the output raster given as <i>glob_rad</i>
+<p>
+The solar irradiation or irradiance raster maps <i>beam_rad</i>, <i>diff_rad</i>,
+<i>refl_rad</i> are computed for a given day <i>day,</i> latitude <i>latin</i>,
+elevation <i>elevin</i>, slope <i>slopein</i> and aspect <i>aspin</i> raster maps.
+For convenience, the output raster given as <i>glob_rad</i>
will output the sum of the three radiation components. The program uses
the Linke atmosphere turbidity factor and ground albedo coefficient.
A default, single value of Linke factor is <i>lin</i>=3.0 and
Modified: grass/trunk/raster/r.sun/rsunlib.c
===================================================================
--- grass/trunk/raster/r.sun/rsunlib.c 2009-08-17 13:19:15 UTC (rev 38771)
+++ grass/trunk/raster/r.sun/rsunlib.c 2009-08-17 13:20:58 UTC (rev 38772)
@@ -183,7 +183,6 @@
costimeAngle = cos(sungeom->timeAngle);
-
lum_Lx = -sungeom->lum_C22 * sin(sungeom->timeAngle);
lum_Ly = sungeom->lum_C11 * costimeAngle + sungeom->lum_C13;
sunVarGeom->sinSolarAltitude =
@@ -240,7 +239,7 @@
inputAngle = (inputAngle >= pi2) ? inputAngle - pi2 : inputAngle;
- delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
+ delt_lat = -0.0001 * cos(inputAngle); /* Arbitrary small distance in latitude */
delt_lon = 0.0001 * sin(inputAngle) / cos(latitude);
newLatitude = (latitude + delt_lat) * rad2deg;
@@ -259,7 +258,6 @@
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;
More information about the grass-commit
mailing list