[GRASS-SVN] r73520 - grass/trunk/imagery/i.eb.hsebal01
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Oct 11 02:24:22 PDT 2018
Author: ychemin
Date: 2018-10-11 02:24:22 -0700 (Thu, 11 Oct 2018)
New Revision: 73520
Modified:
grass/trunk/imagery/i.eb.hsebal01/main.c
Log:
faster/cleaner (all messages appear when verbose level above std)
Modified: grass/trunk/imagery/i.eb.hsebal01/main.c
===================================================================
--- grass/trunk/imagery/i.eb.hsebal01/main.c 2018-10-11 09:22:33 UTC (rev 73519)
+++ grass/trunk/imagery/i.eb.hsebal01/main.c 2018-10-11 09:24:22 UTC (rev 73520)
@@ -43,6 +43,8 @@
return m;
}
+
+
int main(int argc, char *argv[])
{
struct Cell_head cellhd;
@@ -85,8 +87,6 @@
int rowDry = 0, colDry = 0, rowWet = 0, colWet = 0;
/********************************/
-
- /********************************/
xp = yp;
yp = xp;
xmin = ymin;
@@ -135,6 +135,11 @@
input_t0dem->description =
_("Name of altitude corrected surface temperature raster map [K]");
+ input_eact = G_define_standard_option(G_OPT_R_INPUT);
+ input_eact->key = "vapourpressureactual";
+ input_eact->description =
+ _("Name of the actual vapour pressure (e_act) map [KPa]");
+
input_ustar = G_define_option();
input_ustar->key = "frictionvelocitystar";
input_ustar->type = TYPE_DOUBLE;
@@ -145,12 +150,6 @@
_("Value of the height independent friction velocity (u*) [m/s]");
input_ustar->guisection = _("Parameters");
- input_eact = G_define_standard_option(G_OPT_R_INPUT);
- input_eact->key = "vapourpressureactual";
- input_eact->required = YES;
- input_eact->description =
- _("Name of the actual vapour pressure (e_act) map [KPa]");
-
input_row_wet = G_define_option();
input_row_wet->key = "row_wet_pixel";
input_row_wet->type = TYPE_DOUBLE;
@@ -209,7 +208,7 @@
/*If automatic flag, just forget the rest of options */
if (flag2->answer)
- G_message(_("Automatic mode selected"));
+ G_verbose_message(_("Automatic mode selected"));
/*If not automatic & all pixels locations in col/row given */
else if (!flag2->answer &&
input_row_wet->answer &&
@@ -222,9 +221,9 @@
m_col_dry = atof(input_col_dry->answer);
/*If pixels locations are in projected coordinates */
if (flag3->answer)
- G_message(_("Manual wet/dry pixels in image coordinates"));
- G_message(_("Wet Pixel=> x:%f y:%f"), m_col_wet, m_row_wet);
- G_message(_("Dry Pixel=> x:%f y:%f"), m_col_dry, m_row_dry);
+ G_verbose_message(_("Manual wet/dry pixels in image coordinates"));
+ G_verbose_message(_("Wet Pixel=> x:%f y:%f"), m_col_wet, m_row_wet);
+ G_verbose_message(_("Dry Pixel=> x:%f y:%f"), m_col_dry, m_row_dry);
}
/*If not automatic & missing any of the pixel location, Fatal Error */
else {
@@ -278,17 +277,29 @@
/***************************************************/
/* Allocate memory for temporary images */
+
+ /***************************************************/
double **d_Roh, **d_Rah;
+ double **d_z0m, **d_t0dem, **d_eact;
if ((d_Roh = G_alloc_matrix(nrows, ncols)) == NULL)
- G_message("Unable to allocate memory for temporary d_Roh image");
+ G_verbose_message("Unable to allocate memory for temporary d_Roh image");
if ((d_Rah = G_alloc_matrix(nrows, ncols)) == NULL)
- G_message("Unable to allocate memory for temporary d_Rah image");
+ G_verbose_message("Unable to allocate memory for temporary d_Rah image");
+ if ((d_z0m = G_alloc_matrix(nrows, ncols)) == NULL)
+ G_verbose_message("Unable to allocate memory for temporary d_z0m image");
+ if ((d_t0dem = G_alloc_matrix(nrows, ncols)) == NULL)
+ G_verbose_message("Unable to allocate memory for temporary d_t0dem image");
+ if ((d_eact = G_alloc_matrix(nrows, ncols)) == NULL)
+ G_verbose_message("Unable to allocate memory for temporary d_eact image");
/***************************************************/
-
/* MANUAL T0DEM WET/DRY PIXELS */
+ DCELL d_u5 = 0.0;
+ DCELL d_roh1 = 0.0;
+ DCELL d_rah1 = 0.0;
DCELL d_Rn_dry = 0.0, d_g0_dry = 0.0;
+ DCELL d_Rah_dry = 0.0, d_Roh_dry = 0.0;
DCELL d_t0dem_dry = 0.0, d_t0dem_wet = 0.0;
if (flag2->answer) {
@@ -300,74 +311,83 @@
/*********************/
for (row = 0; row < nrows; row++) {
- DCELL d_t0dem;
- DCELL d_z0m;
- DCELL d_eact;
-
G_percent(row, nrows, 2);
Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
Rast_get_d_row(infd_z0m, inrast_z0m, row);
- Rast_get_d_row(infd_eact, inrast_eact, row);
Rast_get_d_row(infd_Rn, inrast_Rn, row);
Rast_get_d_row(infd_g0, inrast_g0, row);
/*process the data */
for (col = 0; col < ncols; col++) {
- d_t0dem = ((DCELL *) inrast_t0dem)[col];
- d_z0m = ((DCELL *) inrast_z0m)[col];
- d_eact = ((DCELL *) inrast_eact)[col];
+ d_t0dem[row][col] = (double) ((DCELL *) inrast_t0dem)[col];
+ d_z0m[row][col] = (double) ((DCELL *) inrast_z0m)[col];
d_Rn = ((DCELL *) inrast_Rn)[col];
d_g0 = ((DCELL *) inrast_g0)[col];
- if (Rast_is_d_null_value(&d_t0dem) ||
- Rast_is_d_null_value(&d_eact) ||
- Rast_is_d_null_value(&d_z0m) ||
- Rast_is_d_null_value(&d_Rn) ||
+ if (Rast_is_d_null_value(&d_Rn) ||
Rast_is_d_null_value(&d_g0)) {
+ d_Roh[row][col] = -999.9;
+ d_Rah[row][col] = -999.9;
/* do nothing */
}
else {
- if (d_t0dem <= 250.0) {
+ if (d_t0dem[row][col] <= 250.0 || d_z0m[row][col] < 0.01) {
+ d_Roh[row][col] = -999.9;
+ d_Rah[row][col] = -999.9;
/* do nothing */
}
else {
d_h0 = d_Rn - d_g0;
- if (d_t0dem < t0dem_min &&
+ d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+ d_rah1 =(1/(d_u5*pow(0.41,2)))*log(5/d_z0m[row][col])*log(5/(d_z0m[row][col]*0.1));
+ d_roh1 =((998-d_eact[row][col])/(d_t0dem[row][col]*2.87))+(d_eact[row][col]/(d_t0dem[row][col]*4.61));
+ if (d_roh1 > 5)
+ d_roh1 = 1.0;
+ else
+ d_roh1 =((1000 - 4.65) / (d_t0dem[row][col] * 2.87))+(4.65 / (d_t0dem[row][col] * 4.61));
+
+ d_Roh[row][col] = d_roh1;
+ d_Rah[row][col] = d_rah1;
+
+ if (d_t0dem[row][col] < t0dem_min &&
d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 0.0 &&
- d_h0 < 100.0) {
- t0dem_min = d_t0dem;
- d_t0dem_wet = d_t0dem;
+ d_h0 < 100.0 && d_roh1 > 0.001 && d_rah1 > 0.001) {
+ t0dem_min = d_t0dem[row][col];
+ d_t0dem_wet = d_t0dem[row][col];
d_Rn_wet = d_Rn;
d_g0_wet = d_g0;
- m_col_wet = col;
- m_row_wet = row;
+ colWet = col;
+ rowWet = row;
}
- if (d_t0dem > t0dem_max &&
+ if (d_t0dem[row][col] > t0dem_max &&
d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 100.0 &&
- d_h0 < 500.0) {
- t0dem_max = d_t0dem;
- d_t0dem_dry = d_t0dem;
+ d_h0 < 500.0 && d_roh1 > 0.001 && d_rah1 > 0.001) {
+ t0dem_max = d_t0dem[row][col];
+ d_t0dem_dry = d_t0dem[row][col];
d_Rn_dry = d_Rn;
d_g0_dry = d_g0;
- m_col_dry = col;
- m_row_dry = row;
+ colDry = col;
+ rowDry = row;
+ d_Roh_dry = d_roh1;
+ d_Rah_dry = d_rah1;
}
}
}
}
}
- G_message("row_wet=%d\tcol_wet=%d", row_wet, col_wet);
- G_message("row_dry=%d\tcol_dry=%d", row_dry, col_dry);
- G_message("g0_wet=%f", d_g0_wet);
- G_message("Rn_wet=%f", d_Rn_wet);
- G_message("LE_wet=%f", d_Rn_wet - d_g0_wet);
- G_message("t0dem_dry=%f", d_t0dem_dry);
- G_message("rnet_dry=%f", d_Rn_dry);
- G_message("g0_dry=%f", d_g0_dry);
- G_message("h0_dry=%f", d_Rn_dry - d_g0_dry);
- G_message("auto config completed");
+ G_verbose_message("row_wet=%d\tcol_wet=%d", rowWet, colWet);
+ G_verbose_message("row_dry=%d\tcol_dry=%d", rowDry, colDry);
+ G_verbose_message("g0_wet=%f", d_g0_wet);
+ G_verbose_message("Rn_wet=%f", d_Rn_wet);
+ G_verbose_message("LE_wet=%f", d_Rn_wet - d_g0_wet);
+ G_verbose_message("t0dem_wet=%f", d_t0dem_wet);
+ G_verbose_message("t0dem_dry=%f", d_t0dem_dry);
+ G_verbose_message("rnet_dry=%f", d_Rn_dry);
+ G_verbose_message("g0_dry=%f", d_g0_dry);
+ G_verbose_message("h0_dry=%f", d_Rn_dry - d_g0_dry);
+ G_verbose_message("Rah_dry=%f", d_Rah_dry);
+ G_verbose_message("Roh_dry=%f", d_Roh_dry);
+ G_verbose_message("auto config completed");
} /* END OF FLAG2 */
- G_message("Passed here");
-
/* MANUAL T0DEM WET/DRY PIXELS */
/*DRY PIXEL */
if (flag3->answer) {
@@ -374,96 +394,97 @@
/*Calculate coordinates of row/col from projected ones */
row = (int)((ymax - m_row_dry) / (double)stepy);
col = (int)((m_col_dry - xmin) / (double)stepx);
- G_message("Dry Pixel | row:%i col:%i", row, col);
+ G_verbose_message("Manual: Dry Pixel | row:%i col:%i", row, col);
+ rowDry = row;
+ colDry = col;
}
else {
- row = (int)m_row_dry;
- col = (int)m_col_dry;
- G_message("Dry Pixel | row:%i col:%i", row, col);
+ row = rowDry;
+ col = colDry;
+ G_verbose_message("Dry Pixel | row:%i col:%i", row, col);
}
- rowDry = row;
- colDry = col;
- Rast_get_d_row(infd_Rn, inrast_Rn, row);
- Rast_get_d_row(infd_g0, inrast_g0, row);
- Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
- d_Rn_dry = ((DCELL *) inrast_Rn)[col];
- d_g0_dry = ((DCELL *) inrast_g0)[col];
- d_t0dem_dry = ((DCELL *) inrast_t0dem)[col];
/*WET PIXEL */
if (flag3->answer) {
/*Calculate coordinates of row/col from projected ones */
row = (int)((ymax - m_row_wet) / (double)stepy);
col = (int)((m_col_wet - xmin) / (double)stepx);
- G_message("Wet Pixel | row:%i col:%i", row, col);
+ G_verbose_message("Manual: Wet Pixel | row:%i col:%i", row, col);
+ rowWet = row;
+ colWet = col;
}
else {
- row = m_row_wet;
- col = m_col_wet;
- G_message("Wet Pixel | row:%i col:%i", row, col);
+ row = rowWet;
+ col = colWet;
+ G_verbose_message("Wet Pixel | row:%i col:%i", row, col);
}
- rowWet = row;
- colWet = col;
- Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
- d_t0dem_wet = ((DCELL *) inrast_t0dem)[col];
/* END OF MANUAL WET/DRY PIXELS */
+
+ /* Extract end-members */
+ Rast_get_d_row(infd_Rn, inrast_Rn, rowDry);
+ Rast_get_d_row(infd_g0, inrast_g0, rowDry);
+ Rast_get_d_row(infd_t0dem, inrast_t0dem, rowDry);
+ d_Rn_dry = ((DCELL *) inrast_Rn)[colDry];
+ d_g0_dry = ((DCELL *) inrast_g0)[colDry];
+ d_t0dem_dry = ((DCELL *) inrast_t0dem)[colDry];
+
+ Rast_get_d_row(infd_t0dem, inrast_t0dem, rowWet);
+ d_t0dem_wet = ((DCELL *) inrast_t0dem)[colWet];
+
double h_dry;
-
h_dry = d_Rn_dry - d_g0_dry;
- G_message("h_dry = %f", h_dry);
- G_message("t0dem_dry = %f", d_t0dem_dry);
- G_message("t0dem_wet = %f", d_t0dem_wet);
+ G_verbose_message("h_dry = %f", h_dry);
+ G_verbose_message("t0dem_dry = %f", d_t0dem_dry);
+ G_verbose_message("t0dem_wet = %f", d_t0dem_wet);
- DCELL d_rah_dry = 0.0;
- DCELL d_roh_dry = 0.0;
+ DCELL d_rah_dry = d_Rah_dry;
+ DCELL d_roh_dry = d_Roh_dry;
- DCELL d_t0dem, d_z0m, d_eact;
- DCELL d_u5;
- DCELL d_roh1;
DCELL d_h1, d_h2, d_h3;
- DCELL d_rah1, d_rah2, d_rah3;
+ DCELL d_rah2, d_rah3;
DCELL d_L, d_x, d_psih, d_psim;
/* INITIALIZATION */
+ /******************************************************/
+ /*If d_rah_dry and d_roh_dry are not auto found, then:*/
+ if(d_Rah_dry==0.0 && d_Roh_dry==0.0){
+ /***************************************************/
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
- /* read a line input maps into buffers */
- Rast_get_d_row(infd_z0m, inrast_z0m, row);
- Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
- Rast_get_d_row(infd_eact, inrast_eact, row);
/* read every cell in the line buffers */
for (col = 0; col < ncols; col++) {
- d_z0m = ((DCELL *) inrast_z0m)[col];
- d_t0dem = ((DCELL *) inrast_t0dem)[col];
- d_eact = ((DCELL *) inrast_eact)[col];
- if (Rast_is_d_null_value(&d_t0dem) ||
- Rast_is_d_null_value(&d_eact) ||
- Rast_is_d_null_value(&d_z0m)) {
+ d_eact[row][col] = (double) ((DCELL *) inrast_eact)[col];
+ d_z0m[row][col] = (double) ((DCELL *) inrast_z0m)[col];
+ if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+ Rast_is_d_null_value(&d_eact[row][col]) ||
+ Rast_is_d_null_value(&d_z0m[row][col])) {
Rast_set_d_null_value(&outrast[col], 1);
d_Roh[row][col] = -999.9;
d_Rah[row][col] = -999.9;
if (row == rowDry && col == colDry) { /*collect dry pix info */
- d_rah_dry = d_rah1;
- d_roh_dry = d_roh1;
- G_message("d_rah_dry=%f d_roh_dry=%f",d_rah_dry,d_roh_dry);
+ d_rah_dry = d_Rah[row][col];
+ d_roh_dry = d_Roh[row][col];
+ G_verbose_message("Init: d_rah_dry=%f d_roh_dry=%f",d_rah_dry,d_roh_dry);
}
}
else {
- d_u5 = (ustar / 0.41) * log(5 / d_z0m);
- d_rah1 =
- (1 / (d_u5 * pow(0.41, 2))) * log(5 / d_z0m) * log(5/(d_z0m*0.1));
- d_roh1 =
- ((998 - d_eact) / (d_t0dem * 2.87)) + (d_eact / (d_t0dem * 4.61));
+ d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+ d_rah1 =(1/(d_u5*pow(0.41,2)))*log(5/d_z0m[row][col])*log(5/(d_z0m[row][col]*0.1));
+ d_roh1 =((998-d_eact[row][col])/(d_t0dem[row][col]*2.87))+(d_eact[row][col]/(d_t0dem[row][col]*4.61));
if (d_roh1 > 5)
d_roh1 = 1.0;
else
- d_roh1 =
- ((1000 - 4.65) / (d_t0dem * 2.87)) +
- (4.65 / (d_t0dem * 4.61));
- if (row == rowDry && col == colDry) { /*collect dry pix info */
+ d_roh1 =((1000 - 4.65) / (d_t0dem[row][col] * 2.87))+(4.65 / (d_t0dem[row][col] * 4.61));
+ if (row == rowDry && col == colDry) {
+ /*collect dry pix info */
d_rah_dry = d_rah1;
d_roh_dry = d_roh1;
- G_message("d_rah_dry=%f d_roh_dry=%f", d_rah_dry,
- d_roh_dry);
+ G_verbose_message("row=%d col=%d", row, col);
+ G_verbose_message("ustar=%f", ustar);
+ G_verbose_message("d_u5=%f", d_u5);
+ G_verbose_message("d_t0dem_dry=%f", d_t0dem[row][col]);
+ G_verbose_message("d_z0m_dry=%f", d_z0m[row][col]);
+ G_verbose_message("d_rah_dry=%f", d_rah_dry);
+ G_verbose_message("d_roh_dry=%f", d_roh_dry);
}
d_Roh[row][col] = d_roh1;
d_Rah[row][col] = d_rah1;
@@ -470,6 +491,7 @@
}
}
}
+ } /*END OF if !d_Rah_dry and !d_Roh_dry*/
DCELL d_dT_dry;
/*Calculate dT_dry */
@@ -486,27 +508,28 @@
a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
b = (sumy - (a * sumx)) / 2.0;
- G_message("d_dT_dry=%f", d_dT_dry);
- G_message("dT1=%f * t0dem + (%f)", a, b);
+ G_verbose_message("d_dT_dry=%f", d_dT_dry);
+ G_verbose_message("dT1=%f * t0dem + (%f)", a, b);
DCELL d_h_dry = 0.0;
if(isnan(a) || isnan(b)){
- G_fatal_error("Delta T Convergence failed, exiting prematurily, please check output");
+ G_free(outrast);
+ Rast_close(outfd);
+ G_fatal_error("Delta T Convergence failed, exiting prematurily, please check output");
}
/* ITERATION 1 */
+ /***************************************************/
+ /***************************************************/
+ //outfd = Rast_open_old(h0, "");
+ //Rast_get_cellhd(h0, "", &cellhd);
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
- /* read a line input maps into buffers */
- Rast_get_d_row(infd_z0m, inrast_z0m, row);
- Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
/* read every cell in the line buffers */
for (col = 0; col < ncols; col++) {
- d_z0m = ((DCELL *) inrast_z0m)[col];
- d_t0dem = ((DCELL *) inrast_t0dem)[col];
d_rah1 = d_Rah[row][col];
d_roh1 = d_Roh[row][col];
- if (Rast_is_d_null_value(&d_t0dem) ||
- Rast_is_d_null_value(&d_z0m)) {
+ if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+ Rast_is_d_null_value(&d_z0m[row][col])) {
Rast_set_d_null_value(&outrast[col], 1);
}
else {
@@ -513,43 +536,55 @@
if (d_rah1 < 1.0)
d_h1 = 0.0;
else
- d_h1 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah1;
+ d_h1 = (1004 * d_roh1) * (a * d_t0dem[row][col] + b) / d_rah1;
if (d_h1 < 0 && d_h1 > -50) {
d_h1 = 0.0;
}
if (d_h1 < -50 || d_h1 > 1000) {
Rast_set_d_null_value(&outrast[col], 1);
+ }else{
+ outrast[col] = d_h1;
+ d_L =-1004*d_roh1*pow(ustar,3)*d_t0dem[row][col]/(d_h1*9.81*0.41);
+ d_x = pow((1 - 16 * (5 / d_L)), 0.25);
+ d_psim =2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
+ 2 * atan(d_x) + 0.5 * M_PI;
+ d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
+ d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+ d_rah2 =(1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m[row][col]) - d_psim)
+ * log((5 / (d_z0m[row][col] * 0.1)) - d_psih);
+ if (row == rowDry && col == colDry) { /*collect dry pix info */
+ d_h1 = (1004 * d_roh1) * (a * d_t0dem[row][col] + b) / d_rah_dry;
+ d_L =-1004*d_roh1*pow(ustar,3)*d_t0dem[row][col]/(d_h1*9.81*0.41);
+ d_x = pow((1 - 16 * (5 / d_L)), 0.25);
+ d_psim =2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
+ 2 * atan(d_x) + 0.5 * M_PI;
+ d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
+ d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+ d_rah2 =(1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m[row][col]) - d_psim)
+ * log((5 / (d_z0m[row][col] * 0.1)) - d_psih);
+ d_rah_dry = d_rah2;
+ d_h_dry = d_h1;
+ G_verbose_message("d_z0m (dry)=%f",d_z0m[row][col]);
+ G_verbose_message("d_rah1 (dry)=%f",d_rah1);
+ G_verbose_message("d_rah2 (dry)=%f",d_rah2);
+ G_verbose_message("d_h1 (dry)=%f",d_h1);
+ }
+ d_Rah[row][col] = d_rah1;
}
- outrast[col] = d_h1;
- d_L =
- -1004 * d_roh1 * pow(ustar,
- 3) * d_t0dem / (d_h1 * 9.81 * 0.41);
- d_x = pow((1 - 16 * (5 / d_L)), 0.25);
- d_psim =
- 2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
- 2 * atan(d_x) + 0.5 * M_PI;
- d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
- d_u5 = (ustar / 0.41) * log(5 / d_z0m);
- d_rah2 =
- (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m) - d_psim)
- * log((5 / (d_z0m * 0.1)) - d_psih);
- if (row == rowDry && col == colDry) { /*collect dry pix info */
- d_rah_dry = d_rah2;
- d_h_dry = d_h1;
- }
- d_Rah[row][col] = d_rah1;
}
}
- Rast_put_d_row(outfd, outrast);
}
/*Calculate dT_dry */
+ G_verbose_message("d_h_dry=%f",d_h_dry);
+ G_verbose_message("d_rah_dry=%f",d_rah_dry);
+ G_verbose_message("d_roh_dry=%f",d_roh_dry);
d_dT_dry = (d_h_dry * d_rah_dry) / (1004 * d_roh_dry);
/*Calculate coefficients for next dT equation */
/* a = (d_dT_dry-0)/(d_t0dem_dry-d_t0dem_wet); */
/* b = (-1.0) * ( a * d_t0dem_wet ); */
- /* G_message("d_dT_dry=%f",d_dT_dry); */
- /* G_message("dT2=%f * t0dem + (%f)", a, b); */
+ /* G_verbose_message("d_dT_dry=%f",d_dT_dry); */
+ /* G_verbose_message("dT2=%f * t0dem + (%f)", a, b); */
sumx = d_t0dem_wet + d_t0dem_dry;
sumy = d_dT_dry + 0.0;
sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
@@ -556,8 +591,8 @@
sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
b = (sumy - (a * sumx)) / 2.0;
- G_message("d_dT_dry=%f", d_dT_dry);
- G_message("dT1=%f * t0dem + (%f)", a, b);
+ G_verbose_message("d_dT_dry=%f", d_dT_dry);
+ G_verbose_message("dT2=%f * t0dem + (%f)", a, b);
if(isnan(a) || isnan(b)){
G_free(outrast);
Rast_close(outfd);
@@ -565,23 +600,16 @@
}
/* ITERATION 2 */
-
/***************************************************/
-
/***************************************************/
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
- /* read a line input maps into buffers */
- Rast_get_d_row(infd_z0m, inrast_z0m, row);
- Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
/* read every cell in the line buffers */
for (col = 0; col < ncols; col++) {
- d_z0m = ((DCELL *) inrast_z0m)[col];
- d_t0dem = ((DCELL *) inrast_t0dem)[col];
d_rah2 = d_Rah[row][col];
d_roh1 = d_Roh[row][col];
- if (Rast_is_d_null_value(&d_t0dem) ||
- Rast_is_d_null_value(&d_z0m)) {
+ if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+ Rast_is_d_null_value(&d_z0m[row][col])) {
Rast_set_d_null_value(&outrast[col], 1);
}
else {
@@ -589,7 +617,7 @@
d_h2 = 0.0;
}
else {
- d_h2 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah2;
+ d_h2 = (1004 * d_roh1) * (a * d_t0dem[row][col] + b) / d_rah2;
}
if (d_h2 < 0 && d_h2 > -50) {
d_h2 = 0.0;
@@ -600,15 +628,15 @@
outrast[col] = d_h2;
d_L =
-1004 * d_roh1 * pow(ustar,
- 3) * d_t0dem / (d_h2 * 9.81 * 0.41);
+ 3) * d_t0dem[row][col] / (d_h2 * 9.81 * 0.41);
d_x = pow((1 - 16 * (5 / d_L)), 0.25);
d_psim = 2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
2 * atan(d_x) + 0.5 * M_PI;
d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
- d_u5 = (ustar / 0.41) * log(5 / d_z0m);
+ d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
d_rah3 =
- (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m) -
- d_psim) * log((5 / (d_z0m * 0.1)) - d_psih);
+ (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m[row][col]) -
+ d_psim) * log((5 / (d_z0m[row][col] * 0.1)) - d_psih);
if (row == rowDry && col == colDry) { /*collect dry pix info */
d_rah_dry = d_rah3;
d_h_dry = d_h2;
@@ -616,7 +644,6 @@
d_Rah[row][col] = d_rah2;
}
}
- Rast_put_d_row(outfd, outrast);
}
/*Calculate dT_dry */
@@ -624,8 +651,8 @@
/*Calculate coefficients for next dT equation */
/* a = (d_dT_dry-0)/(d_t0dem_dry-d_t0dem_wet); */
/* b = (-1.0) * ( a * d_t0dem_wet ); */
- /* G_message("d_dT_dry=%f",d_dT_dry); */
- /* G_message("dT3=%f * t0dem + (%f)", a, b); */
+ /* G_verbose_message("d_dT_dry=%f",d_dT_dry); */
+ /* G_verbose_message("dT3=%f * t0dem + (%f)", a, b); */
sumx = d_t0dem_wet + d_t0dem_dry;
sumy = d_dT_dry + 0.0;
sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
@@ -632,8 +659,8 @@
sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
b = (sumy - (a * sumx)) / 2.0;
- G_message("d_dT_dry=%f", d_dT_dry);
- G_message("dT1=%f * t0dem + (%f)", a, b);
+ G_verbose_message("d_dT_dry=%f", d_dT_dry);
+ G_verbose_message("dT3=%f * t0dem + (%f)", a, b);
if(isnan(a) || isnan(b)){
G_free(outrast);
Rast_close(outfd);
@@ -641,24 +668,15 @@
}
/* ITERATION 3 */
-
/***************************************************/
-
/***************************************************/
-
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
- /* read a line input maps into buffers */
- Rast_get_d_row(infd_z0m, inrast_z0m, row);
- Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
/* read every cell in the line buffers */
for (col = 0; col < ncols; col++) {
- d_z0m = ((DCELL *) inrast_z0m)[col];
- d_t0dem = ((DCELL *) inrast_t0dem)[col];
d_rah3 = d_Rah[row][col];
- d_roh1 = d_Roh[row][col];
- if (Rast_is_d_null_value(&d_t0dem) ||
- Rast_is_d_null_value(&d_z0m)) {
+ if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+ Rast_is_d_null_value(&d_z0m[row][col])) {
Rast_set_d_null_value(&outrast[col], 1);
}
else {
@@ -666,7 +684,7 @@
d_h3 = 0.0;
}
else {
- d_h3 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah3;
+ d_h3 = (1004 * d_Roh[row][col]) * (a * d_t0dem[row][col] + b) / d_rah3;
}
if (d_h3 < 0 && d_h3 > -50) {
d_h3 = 0.0;
@@ -673,14 +691,15 @@
}
if (d_h3 < -50 || d_h3 > 1000) {
Rast_set_d_null_value(&outrast[col], 1);
+ } else {
+ outrast[col] = d_h3;
}
- outrast[col] = d_h3;
}
}
Rast_put_d_row(outfd, outrast);
}
-
-
+ G_free(inrast_eact);
+ Rast_close(infd_eact);
G_free(inrast_z0m);
Rast_close(infd_z0m);
G_free(inrast_t0dem);
More information about the grass-commit
mailing list