diff --git a/raster/r.horizon/TODO b/raster/r.horizon/TODO deleted file mode 100644 index e31cd59d614..00000000000 --- a/raster/r.horizon/TODO +++ /dev/null @@ -1,11 +0,0 @@ -TODO - -Probably the sun position calculation should be replaced -with - - G_calc_solar_position() - -currently used in r.sunmask. G_calc_solar_position() is based on -solpos.c from NREL and very precise. - -MN 4/2001 diff --git a/raster/r.horizon/local_proto.h b/raster/r.horizon/local_proto.h deleted file mode 100644 index f2b9232d476..00000000000 --- a/raster/r.horizon/local_proto.h +++ /dev/null @@ -1,23 +0,0 @@ -/* sun11.c */ -int INPUT(void); -int OUTGR(void); -double amax1(double, double); -double amin1(double, double); -int min(int, int); -int max(int, int); -void com_par(void); -double lumcline2(void); -double joules2(void); -int quadrant(void); -double coef_of_line(void); -void new_point_x(int, double *, double *, double *); -void new_point_y(int, double *, double *, double *); -void which_one(double, double, double, double, double, double); -int combine_x(int, int, int, int); -int combine_y(int, int, int, int); -int vertex(int, int); -int mesh_vertex(void); -int mesh_line(void); -void calculate(void); -double com_sol_const(void); -double com_declin(int); diff --git a/raster/r.horizon/main.c b/raster/r.horizon/main.c index 5624dc5db2f..d434f119e26 100644 --- a/raster/r.horizon/main.c +++ b/raster/r.horizon/main.c @@ -12,6 +12,8 @@ This program was written in 2006 by Thomas Huld and Tomas Cebecauer, Joint Research Centre of the European Commission, based on bits of the r.sun module by Jaro Hofierka +Program was refactored by Anna Petrasova to remove most global variables. + *******************************************************************************/ /* * This program is free software; you can redistribute it and/or @@ -39,112 +41,101 @@ module by Jaro Hofierka #include #include -#define WHOLE_RASTER 1 -#define SINGLE_POINT 0 -#define RAD (180. / M_PI) -#define DEG ((M_PI) / 180.) -#define EARTHRADIUS 6371000. -#define UNDEF 0. /* undefined value for terrain aspect */ -#define UNDEFZ -9999. /* undefined value for elevation */ -#define BIG 1.e20 -#define SMALL 1.e-20 -#define EPS 1.e-4 -#define DIST "1.0" -#define DEGREEINMETERS 111120. /* 1852m/nm * 60nm/degree = 111120 m/deg */ -#define TANMINANGLE 0.008727 /* tan of minimum horizon angle (0.5 deg) */ - -#define AMAX1(arg1, arg2) ((arg1) >= (arg2) ? (arg1) : (arg2)) +#define WHOLE_RASTER 1 +#define SINGLE_POINT 0 +#define RAD (180. / M_PI) +#define DEG ((M_PI) / 180.) +#define EARTHRADIUS 6371000. +#define UNDEF 0. /* undefined value for terrain aspect */ +#define UNDEFZ -9999. /* undefined value for elevation */ +#define BIG 1.e20 +#define SMALL 1.e-20 +#define EPS 1.e-4 +#define DIST "1.0" +#define DEGREEINMETERS 111120. /* 1852m/nm * 60nm/degree = 111120 m/deg */ +#define TANMINANGLE 0.008727 /* tan of minimum horizon angle (0.5 deg) */ + #define DISTANCE1(x1, x2, y1, y2) \ (sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))) -FILE *fp; - const double pihalf = M_PI * 0.5; const double twopi = M_PI * 2.; const double invEarth = 1. / EARTHRADIUS; const double deg2rad = M_PI / 180.; const double rad2deg = 180. / M_PI; -const double minAngle = DEG; -const char *elevin; -const char *horizon = NULL; -const char *mapset = NULL; -const char *per; -char *shad_filename; -char *outfile; - -struct Cell_head cellhd; -struct Key_value *in_proj_info, *in_unit_info; struct pj_info iproj, oproj, tproj; - -struct Cell_head new_cellhd; -double bufferZone = 0., ebufferZone = 0., wbufferZone = 0., nbufferZone = 0., - sbufferZone = 0.; - -int INPUT(void); -int OUTGR(int numrows, int numcols); -double amax1(double, double); -double amin1(double, double); -int min(int, int); -int max(int, int); -void com_par(void); -int is_shadow(void); -double horizon_height(void); -void calculate_shadow(void); -double calculate_shadow_onedirection(double shadow_angle); - -int new_point(void); -double searching(void); -int test_low_res(void); - -/*void where_is_point(); - void cube(int, int); - */ - -void calculate(double xcoord, double ycoord, int buffer_e, int buffer_w, - int buffer_s, int buffer_n); - -int ip, jp, ip100, jp100; -int n, m, m100, n100; -int degreeOutput, compassOutput = FALSE; float **z, **z100, **horizon_raster; -double stepx, stepy, stepxhalf, stepyhalf, stepxy, xp, yp, op, dp, xg0, xx0, - yg0, yy0, deltx, delty; -double invstepx, invstepy, distxy; -double offsetx, offsety; -double single_direction; - -/*int arrayNumInt; */ -double xmin, xmax, ymin, ymax, zmax = 0.; -int d, day, tien = 0; - -double length, maxlength = BIG, zmult = 1.0, dist; -double fixedMaxLength = BIG, step = 0.0, start = 0.0, end = 0.0; -char *tt, *lt; -double z_orig, zp; -double h0, tanh0, angle; -double stepsinangle, stepcosangle, sinangle, cosangle, distsinangle, - distcosangle; -double TOLER; -const char *str_step; - -int mode; -int isMode(void) -{ - return mode; -} - -void setMode(int val) -{ - mode = val; -} - -int ll_correction = FALSE; -double coslatsq; +bool ll_correction = false; + +typedef struct { + double xg0, yg0; + double z_orig; + double coslatsq; + double maxlength; +} OriginPoint; + +typedef struct { + double stepsinangle, stepcosangle; + double sinangle, cosangle; + double distsinangle, distcosangle; +} OriginAngle; + +typedef struct { + double xx0, yy0; + int ip, jp; + int ip100, jp100; + double zp; +} SearchPoint; + +typedef struct { + double tanh0; + double length; +} HorizonProperties; + +typedef struct { + int n, m, m100, n100; + double stepx, stepy, stepxy; + double invstepx, invstepy; + double offsetx, offsety; + double distxy; + double xmin, xmax, ymin, ymax, zmax; +} Geometry; + +typedef struct { + int degreeOutput; + int compassOutput; + double fixedMaxLength; + double start, end, step; + double single_direction; + const char *str_step; + const char *horizon_basename; +} Settings; + +int INPUT(Geometry *geometry, const char *elevin); +int OUTGR(const Settings *settings, char *shad_filename, + struct Cell_head *cellhd); +void com_par(const Geometry *geometry, OriginAngle *origin_angle, double angle, + double xp, double yp); +double horizon_height(const Geometry *geometry, const OriginPoint *origin_point, + const OriginAngle *origin_angle); +void calculate_point_mode(const Settings *settings, const Geometry *geometry, + double xcoord, double ycoord, FILE *fp); +int new_point(const Geometry *geometry, const OriginPoint *origin_point, + const OriginAngle *origin_angle, SearchPoint *search_point, + HorizonProperties *horizon); +int test_low_res(const Geometry *geometry, const OriginPoint *origin_point, + const OriginAngle *origin_angle, SearchPoint *search_point, + const HorizonProperties *horizon); +void calculate_raster_mode(const Settings *settings, const Geometry *geometry, + struct Cell_head *cellhd, + struct Cell_head *new_cellhd, int buffer_e, + int buffer_w, int buffer_s, int buffer_n, + double bufferZone); /* why not use G_distance() here which switches to geodesic/great circle distance as needed? */ -double distance(double x1, double x2, double y1, double y2) +double distance(double x1, double x2, double y1, double y2, double coslatsq) { if (ll_correction) { return DEGREEINMETERS * @@ -318,52 +309,53 @@ int main(int argc, char *argv[]) if (G_parser(argc, argv)) exit(EXIT_FAILURE); + struct Cell_head cellhd; + struct Cell_head new_cellhd; G_get_set_window(&cellhd); + Geometry geometry; - stepx = cellhd.ew_res; - stepy = cellhd.ns_res; - stepxhalf = stepx / 2.; - stepyhalf = stepy / 2.; - invstepx = 1. / stepx; - invstepy = 1. / stepy; + geometry.stepx = cellhd.ew_res; + geometry.stepy = cellhd.ns_res; + geometry.invstepx = 1. / geometry.stepx; + geometry.invstepy = 1. / geometry.stepy; /* offsetx = 2. * invstepx; offsety = 2. * invstepy; offsetx = 0.5*stepx; offsety = 0.5*stepy; */ - offsetx = 0.5; - offsety = 0.5; + geometry.offsetx = 0.5; + geometry.offsety = 0.5; - n /*n_cols */ = cellhd.cols; - m /*n_rows */ = cellhd.rows; + geometry.n /*n_cols */ = cellhd.cols; + geometry.m /*n_rows */ = cellhd.rows; - n100 = ceil(n / 100.); - m100 = ceil(m / 100.); + geometry.n100 = ceil(geometry.n / 100.); + geometry.m100 = ceil(geometry.m / 100.); - xmin = cellhd.west; - ymin = cellhd.south; - xmax = cellhd.east; - ymax = cellhd.north; - deltx = fabs(cellhd.east - cellhd.west); - delty = fabs(cellhd.north - cellhd.south); + geometry.xmin = cellhd.west; + geometry.ymin = cellhd.south; + geometry.xmax = cellhd.east; + geometry.ymax = cellhd.north; - degreeOutput = flag.degreeOutput->answer; - compassOutput = flag.compassOutput->answer; + Settings settings; + settings.degreeOutput = flag.degreeOutput->answer; + settings.compassOutput = flag.compassOutput->answer; if (G_projection() == PROJECTION_LL) G_important_message(_("Note: In latitude-longitude coordinate system " "specify buffers in degree unit")); - elevin = parm.elevin->answer; - + const char *elevin = parm.elevin->answer; + FILE *fp = NULL; + int mode; if (parm.coord->answer == NULL) { G_debug(1, "Setting mode: WHOLE_RASTER"); - setMode(WHOLE_RASTER); + mode = WHOLE_RASTER; } else { G_debug(1, "Setting mode: SINGLE_POINT"); - setMode(SINGLE_POINT); + mode = SINGLE_POINT; if (sscanf(parm.coord->answer, "%lf,%lf", &xcoord, &ycoord) != 2) { G_fatal_error(_( "Can't read the coordinates from the \"coordinate\" option.")); @@ -382,7 +374,7 @@ int main(int argc, char *argv[]) */ /* Open ASCII file for output or stdout */ - outfile = parm.output->answer; + char *outfile = parm.output->answer; if ((strcmp("-", outfile)) == 0) { fp = stdout; @@ -390,11 +382,15 @@ int main(int argc, char *argv[]) else if (NULL == (fp = fopen(outfile, "w"))) G_fatal_error(_("Unable to open file <%s>"), outfile); } - + settings.single_direction = 0; if (parm.direction->answer != NULL) - sscanf(parm.direction->answer, "%lf", &single_direction); + sscanf(parm.direction->answer, "%lf", &settings.single_direction); - if (WHOLE_RASTER == isMode()) { + settings.step = 0; + settings.start = 0; + settings.end = 0; + settings.horizon_basename = NULL; + if (WHOLE_RASTER == mode) { if ((parm.direction->answer == NULL) && (parm.step->answer == NULL)) { G_fatal_error(_("You didn't specify a direction value or step " "size. Aborting.")); @@ -404,29 +400,30 @@ int main(int argc, char *argv[]) G_fatal_error( _("You didn't specify a horizon raster name. Aborting.")); } - horizon = parm.horizon->answer; + settings.horizon_basename = parm.horizon->answer; if (parm.step->answer != NULL) { - str_step = parm.step->answer; - sscanf(parm.step->answer, "%lf", &step); + settings.str_step = parm.step->answer; + sscanf(parm.step->answer, "%lf", &settings.step); } else { - step = 0; - str_step = "0"; + settings.step = 0; + settings.str_step = "0"; } - sscanf(parm.start->answer, "%lf", &start); - sscanf(parm.end->answer, "%lf", &end); - if (start < 0.0) { + sscanf(parm.start->answer, "%lf", &settings.start); + sscanf(parm.end->answer, "%lf", &settings.end); + if (settings.start < 0.0) { G_fatal_error( _("Negative values of start angle are not allowed. Aborting.")); } - if (end < 0.0 || end > 360.0) { + if (settings.end < 0.0 || settings.end > 360.0) { G_fatal_error(_("End angle is not between 0 and 360. Aborting.")); } - if (start >= end) { + if (settings.start >= settings.end) { G_fatal_error( _("You specify a start grater than the end angle. Aborting.")); } - G_debug(1, "Angle step: %g, start: %g, end: %g", step, start, end); + G_debug(1, "Angle step: %g, start: %g, end: %g", settings.step, + settings.start, settings.end); } else { @@ -434,13 +431,14 @@ int main(int argc, char *argv[]) G_fatal_error( _("You didn't specify an angle step size. Aborting.")); } - sscanf(parm.step->answer, "%lf", &step); + sscanf(parm.step->answer, "%lf", &settings.step); } - if (step == 0.0) { - step = 360.; + if (settings.step == 0.0) { + settings.step = 360.; } - + double bufferZone = 0., ebufferZone = 0., wbufferZone = 0., + nbufferZone = 0., sbufferZone = 0.; if (parm.bufferzone->answer != NULL) { if (sscanf(parm.bufferzone->answer, "%lf", &bufferZone) != 1) G_fatal_error(_("Could not read bufferzone size. Aborting.")); @@ -470,11 +468,14 @@ int main(int argc, char *argv[]) _("north")); } + settings.fixedMaxLength = BIG; if (parm.maxdistance->answer != NULL) { - if (sscanf(parm.maxdistance->answer, "%lf", &fixedMaxLength) != 1) + if (sscanf(parm.maxdistance->answer, "%lf", &settings.fixedMaxLength) != + 1) G_fatal_error(_("Could not read maximum distance. Aborting.")); } - G_debug(1, "Using maxdistance %f", fixedMaxLength); /* predefined as BIG */ + G_debug(1, "Using maxdistance %f", + settings.fixedMaxLength); /* predefined as BIG */ /* TODO: fixing BIG, there is a bug with distant mountains not being seen: attempt to contrain to current region @@ -482,14 +483,11 @@ int main(int argc, char *argv[]) fixedMaxLength = (fixedMaxLength < AMAX1(deltx, delty)) ? fixedMaxLength : AMAX1(deltx, delty); G_debug(1,"Using maxdistance %f", fixedMaxLength); */ - - sscanf(parm.dist->answer, "%lf", &dist); - if (dist < 0.5 || dist > 1.5) + sscanf(parm.dist->answer, "%lf", &geometry.distxy); + if (geometry.distxy < 0.5 || geometry.distxy > 1.5) G_fatal_error(_("The distance value must be 0.5-1.5. Aborting.")); - stepxy = dist * 0.5 * (stepx + stepy); - distxy = dist; - TOLER = stepxy * EPS; + geometry.stepxy = geometry.distxy * 0.5 * (geometry.stepx + geometry.stepy); if (bufferZone > 0. || ebufferZone > 0. || wbufferZone > 0. || sbufferZone > 0. || nbufferZone > 0.) { @@ -505,31 +503,30 @@ int main(int argc, char *argv[]) nbufferZone = bufferZone; /* adjust buffer to multiples of resolution */ - ebufferZone = (int)(ebufferZone / stepx) * stepx; - wbufferZone = (int)(wbufferZone / stepx) * stepx; - sbufferZone = (int)(sbufferZone / stepy) * stepy; - nbufferZone = (int)(nbufferZone / stepy) * stepy; + ebufferZone = (int)(ebufferZone / geometry.stepx) * geometry.stepx; + wbufferZone = (int)(wbufferZone / geometry.stepx) * geometry.stepx; + sbufferZone = (int)(sbufferZone / geometry.stepy) * geometry.stepy; + nbufferZone = (int)(nbufferZone / geometry.stepy) * geometry.stepy; - new_cellhd.rows += (int)((nbufferZone + sbufferZone) / stepy); - new_cellhd.cols += (int)((ebufferZone + wbufferZone) / stepx); + new_cellhd.rows += (int)((nbufferZone + sbufferZone) / geometry.stepy); + new_cellhd.cols += (int)((ebufferZone + wbufferZone) / geometry.stepx); new_cellhd.north += nbufferZone; new_cellhd.south -= sbufferZone; new_cellhd.east += ebufferZone; new_cellhd.west -= wbufferZone; - xmin = new_cellhd.west; - ymin = new_cellhd.south; - xmax = new_cellhd.east; - ymax = new_cellhd.north; - deltx = fabs(new_cellhd.east - new_cellhd.west); - delty = fabs(new_cellhd.north - new_cellhd.south); + geometry.xmin = new_cellhd.west; + geometry.ymin = new_cellhd.south; + geometry.xmax = new_cellhd.east; + geometry.ymax = new_cellhd.north; - n /* n_cols */ = new_cellhd.cols; - m /* n_rows */ = new_cellhd.rows; - G_debug(1, "%lf %lf %lf %lf \n", ymax, ymin, xmin, xmax); - n100 = ceil(n / 100.); - m100 = ceil(m / 100.); + geometry.n /* n_cols */ = new_cellhd.cols; + geometry.m /* n_rows */ = new_cellhd.rows; + G_debug(1, "%lf %lf %lf %lf \n", geometry.ymax, geometry.ymin, + geometry.xmin, geometry.xmax); + geometry.n100 = ceil(geometry.n / 100.); + geometry.m100 = ceil(geometry.m / 100.); Rast_set_window(&new_cellhd); } @@ -556,44 +553,46 @@ int main(int argc, char *argv[]) /**********end of parser - ******************************/ - INPUT(); - G_debug(1, "calculate() starts..."); - calculate(xcoord, ycoord, (int)(ebufferZone / stepx), - (int)(wbufferZone / stepx), (int)(sbufferZone / stepy), - (int)(nbufferZone / stepy)); + INPUT(&geometry, elevin); + if (mode == SINGLE_POINT) { + /* Calculate the horizon for one single point */ + calculate_point_mode(&settings, &geometry, xcoord, ycoord, fp); + } + else { + calculate_raster_mode(&settings, &geometry, &cellhd, &new_cellhd, + (int)(ebufferZone / geometry.stepx), + (int)(wbufferZone / geometry.stepx), + (int)(sbufferZone / geometry.stepy), + (int)(nbufferZone / geometry.stepy), bufferZone); + } exit(EXIT_SUCCESS); } /**********************end of main.c*****************/ -int INPUT(void) +int INPUT(Geometry *geometry, const char *elevin) { - FCELL *cell1; - int fd1, row, row_rev; - int l, i, j, k; - int lmax, kmax; + FCELL *cell1 = Rast_allocate_f_buf(); - cell1 = Rast_allocate_f_buf(); + z = (float **)G_malloc(sizeof(float *) * (geometry->m)); + z100 = (float **)G_malloc(sizeof(float *) * (geometry->m100)); - z = (float **)G_malloc(sizeof(float *) * (m)); - z100 = (float **)G_malloc(sizeof(float *) * (m100)); - - for (l = 0; l < m; l++) { - z[l] = (float *)G_malloc(sizeof(float) * (n)); + for (int l = 0; l < geometry->m; l++) { + z[l] = (float *)G_malloc(sizeof(float) * (geometry->n)); } - for (l = 0; l < m100; l++) { - z100[l] = (float *)G_malloc(sizeof(float) * (n100)); + for (int l = 0; l < geometry->m100; l++) { + z100[l] = (float *)G_malloc(sizeof(float) * (geometry->n100)); } /*read Z raster */ - fd1 = Rast_open_old(elevin, ""); + int fd1 = Rast_open_old(elevin, ""); - for (row = 0; row < m; row++) { + for (int row = 0; row < geometry->m; row++) { Rast_get_f_row(fd1, cell1, row); - for (j = 0; j < n; j++) { - row_rev = m - row - 1; + for (int j = 0; j < geometry->n; j++) { + int row_rev = geometry->m - row - 1; if (!Rast_is_f_null_value(cell1 + j)) z[row_rev][j] = (float)cell1[j]; @@ -604,46 +603,47 @@ int INPUT(void) Rast_close(fd1); /* create low resolution array 100 */ - for (i = 0; i < m100; i++) { - lmax = (i + 1) * 100; - if (lmax > m) - lmax = m; - - for (j = 0; j < n100; j++) { - zmax = SMALL; - kmax = (j + 1) * 100; - if (kmax > n) - kmax = n; - for (l = (i * 100); l < lmax; l++) { - for (k = (j * 100); k < kmax; k++) { - zmax = amax1(zmax, z[l][k]); + for (int i = 0; i < geometry->m100; i++) { + int lmax = (i + 1) * 100; + if (lmax > geometry->m) + lmax = geometry->m; + + for (int j = 0; j < geometry->n100; j++) { + geometry->zmax = SMALL; + int kmax = (j + 1) * 100; + if (kmax > geometry->n) + kmax = geometry->n; + for (int l = (i * 100); l < lmax; l++) { + for (int k = (j * 100); k < kmax; k++) { + geometry->zmax = MAX(geometry->zmax, z[l][k]); } } - z100[i][j] = zmax; + z100[i][j] = geometry->zmax; G_debug(3, "%d %d %lf\n", i, j, z100[i][j]); } } /* find max Z */ - for (i = 0; i < m; i++) { - for (j = 0; j < n; j++) { - zmax = amax1(zmax, z[i][j]); + for (int i = 0; i < geometry->m; i++) { + for (int j = 0; j < geometry->n; j++) { + geometry->zmax = MAX(geometry->zmax, z[i][j]); } } return 1; } -int OUTGR(int numrows, int numcols) +int OUTGR(const Settings *settings, char *shad_filename, + struct Cell_head *cellhd) { FCELL *cell1 = NULL; - int fd1 = 0; - int i, iarc, j; - Rast_set_window(&cellhd); + int numrows = cellhd->rows; + int numcols = cellhd->cols; + Rast_set_window(cellhd); - if (horizon != NULL) { + if (settings->horizon_basename != NULL) { cell1 = Rast_allocate_f_buf(); fd1 = Rast_open_fp_new(shad_filename); } @@ -656,11 +656,11 @@ int OUTGR(int numrows, int numcols) G_fatal_error(_("OOPS: cols changed from %d to %d"), numcols, Rast_window_cols()); - for (iarc = 0; iarc < numrows; iarc++) { - i = numrows - iarc - 1; + for (int iarc = 0; iarc < numrows; iarc++) { + int i = numrows - iarc - 1; - if (horizon != NULL) { - for (j = 0; j < numcols; j++) { + if (settings->horizon_basename != NULL) { + for (int j = 0; j < numcols; j++) { if (horizon_raster[i][j] == UNDEFZ) Rast_set_f_null_value(cell1 + j, 1); else @@ -675,189 +675,119 @@ int OUTGR(int numrows, int numcols) return 1; } -double amax1(double arg1, double arg2) -{ - double res; - - if (arg1 >= arg2) { - res = arg1; - } - else { - res = arg2; - } - return res; -} +/**********************************************************/ -double amin1(double arg1, double arg2) +void com_par(const Geometry *geometry, OriginAngle *origin_angle, double angle, + double xp, double yp) { - double res; - - if (arg1 <= arg2) { - res = arg1; - } - else { - res = arg2; + double longitude = xp; + double latitude = yp; + if (G_projection() != PROJECTION_LL) { + if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD, &longitude, &latitude, + NULL) < 0) + G_fatal_error(_("Error in %s"), "GPJ_transform()"); } - return res; -} + latitude *= deg2rad; + longitude *= deg2rad; -int min(int arg1, int arg2) -{ - int res; + double delt_lat = + -0.0001 * cos(angle); /* Arbitrary small distance in latitude */ + double delt_lon = 0.0001 * sin(angle) / cos(latitude); - if (arg1 <= arg2) { - res = arg1; - } - else { - res = arg2; - } - return res; -} + latitude = (latitude + delt_lat) * rad2deg; + longitude = (longitude + delt_lon) * rad2deg; -int max(int arg1, int arg2) -{ - int res; - - if (arg1 >= arg2) { - res = arg1; - } - else { - res = arg2; + if (G_projection() != PROJECTION_LL) { + if (GPJ_transform(&iproj, &oproj, &tproj, PJ_INV, &longitude, &latitude, + NULL) < 0) + G_fatal_error(_("Error in %s"), "GPJ_transform()"); } - return res; -} + double delt_east = longitude - xp; + double delt_nor = latitude - yp; -/**********************************************************/ + double delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor); -void com_par(void) -{ - if (fabs(sinangle) < 0.0000001) { - sinangle = 0.; + origin_angle->sinangle = delt_nor / delt_dist; + origin_angle->cosangle = delt_east / delt_dist; + + if (fabs(origin_angle->sinangle) < 0.0000001) { + origin_angle->sinangle = 0.; } - if (fabs(cosangle) < 0.0000001) { - cosangle = 0.; + if (fabs(origin_angle->cosangle) < 0.0000001) { + origin_angle->cosangle = 0.; } - distsinangle = 32000; - distcosangle = 32000; + origin_angle->distsinangle = 32000; + origin_angle->distcosangle = 32000; - if (sinangle != 0.) { - distsinangle = 100. / (distxy * sinangle); + if (origin_angle->sinangle != 0.) { + origin_angle->distsinangle = + 100. / (geometry->distxy * origin_angle->sinangle); } - if (cosangle != 0.) { - distcosangle = 100. / (distxy * cosangle); + if (origin_angle->cosangle != 0.) { + origin_angle->distcosangle = + 100. / (geometry->distxy * origin_angle->cosangle); } - stepsinangle = stepxy * sinangle; - stepcosangle = stepxy * cosangle; -} - -double horizon_height(void) -{ - double height; - - tanh0 = 0.; - length = 0; - zp = z_orig; - - height = searching(); - - xx0 = xg0; - yy0 = yg0; - - return height; + origin_angle->stepsinangle = geometry->stepxy * origin_angle->sinangle; + origin_angle->stepcosangle = geometry->stepxy * origin_angle->cosangle; } -double calculate_shadow_onedirection(double shadow_angle) +void calculate_point_mode(const Settings *settings, const Geometry *geometry, + double xcoord, double ycoord, FILE *fp) { - shadow_angle = horizon_height(); - - return shadow_angle; -} - -void calculate_shadow(void) -{ - double dfr_rad; - - int i; - int printCount; - double shadow_angle; - double printangle; - double sx, sy; - double xp, yp; - double latitude, longitude; - double delt_lat, delt_lon; - double delt_east, delt_nor; - double delt_dist; + /* + xg0 = xx0 = (double)xcoord * stepx; + yg0 = yy0 = (double)ycoord * stepy; + xg0 = xx0 = xcoord -0.5*stepx -xmin; + yg0 = yy0 = ycoord -0.5*stepy-ymin; + xg0 = xx0 = xindex*stepx -0.5*stepx; + yg0 = yy0 = yindex*stepy -0.5*stepy; + */ + OriginPoint origin_point; + int xindex = (int)((xcoord - geometry->xmin) / geometry->stepx); + int yindex = (int)((ycoord - geometry->ymin) / geometry->stepy); + origin_point.xg0 = xindex * geometry->stepx; + origin_point.yg0 = yindex * geometry->stepy; + origin_point.coslatsq = 0; + if ((G_projection() == PROJECTION_LL)) { + ll_correction = true; + } + if (ll_correction) { + double coslat = cos(deg2rad * (geometry->ymin + origin_point.yg0)); + origin_point.coslatsq = coslat * coslat; + } - double angle; + origin_point.z_orig = z[yindex][xindex]; + G_debug(1, "yindex: %d, xindex %d, z_orig %.2f", yindex, xindex, + origin_point.z_orig); - printCount = 360. / fabs(step); + int printCount = 360. / fabs(settings->step); if (printCount < 1) printCount = 1; + double dfr_rad = settings->step * deg2rad; - dfr_rad = step * deg2rad; - - xp = xmin + xx0; - yp = ymin + yy0; + double xp = geometry->xmin + origin_point.xg0; + double yp = geometry->ymin + origin_point.yg0; - angle = (single_direction * deg2rad) + pihalf; - printangle = single_direction; + double angle = (settings->single_direction * deg2rad) + pihalf; + double printangle = settings->single_direction; - maxlength = fixedMaxLength; + origin_point.maxlength = settings->fixedMaxLength; fprintf(fp, "azimuth,horizon_height\n"); - for (i = 0; i < printCount; i++) { - - ip = jp = 0; - - sx = xx0 * invstepx; - sy = yy0 * invstepy; - ip100 = floor(sx / 100.); - jp100 = floor(sy / 100.); - - if ((G_projection() != PROJECTION_LL)) { - - longitude = xp; - latitude = yp; - - if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD, &longitude, - &latitude, NULL) < 0) - G_fatal_error(_("Error in %s"), "GPJ_transform()"); - } - else { /* ll projection */ - latitude = yp; - longitude = xp; - } - latitude *= deg2rad; - longitude *= deg2rad; - - delt_lat = -0.0001 * cos(angle); - delt_lon = 0.0001 * sin(angle) / cos(latitude); - - latitude = (latitude + delt_lat) * rad2deg; - longitude = (longitude + delt_lon) * rad2deg; - - if (GPJ_transform(&iproj, &oproj, &tproj, PJ_INV, &longitude, &latitude, - NULL) < 0) - G_fatal_error(_("Error in %s"), "GPJ_transform()"); - - delt_east = longitude - xp; - delt_nor = latitude - yp; - - delt_dist = sqrt(delt_east * delt_east + delt_nor * delt_nor); + for (int i = 0; i < printCount; i++) { + OriginAngle origin_angle; + com_par(geometry, &origin_angle, angle, xp, yp); - sinangle = delt_nor / delt_dist; - cosangle = delt_east / delt_dist; - com_par(); + double shadow_angle = + horizon_height(geometry, &origin_point, &origin_angle); - shadow_angle = horizon_height(); - - if (degreeOutput) { + if (settings->degreeOutput) { shadow_angle *= rad2deg; } - if (compassOutput) { + if (settings->compassOutput) { double tmpangle; tmpangle = 360. - printangle + 90.; @@ -870,7 +800,7 @@ void calculate_shadow(void) } angle += dfr_rad; - printangle += step; + printangle += settings->step; if (angle < 0.) angle += twopi; @@ -882,44 +812,45 @@ void calculate_shadow(void) else if (printangle > 360.) printangle -= 360; } /* end of for loop over angles */ + fclose(fp); } /*////////////////////////////////////////////////////////////////////// */ -int new_point(void) +int new_point(const Geometry *geometry, const OriginPoint *origin_point, + const OriginAngle *origin_angle, SearchPoint *search_point, + HorizonProperties *horizon) { - int iold, jold; - int succes = 1, succes2 = 1; - double sx, sy; - double dx, dy; - - iold = ip; - jold = jp; + int iold = search_point->ip; + int jold = search_point->jp; - while (succes) { - yy0 += stepsinangle; - xx0 += stepcosangle; + while (TRUE) { + search_point->yy0 += origin_angle->stepsinangle; + search_point->xx0 += origin_angle->stepcosangle; /* offset 0.5 cell size to get the right cell i, j */ - sx = xx0 * invstepx + offsetx; - sy = yy0 * invstepy + offsety; - ip = (int)sx; - jp = (int)sy; + double sx = search_point->xx0 * geometry->invstepx + geometry->offsetx; + double sy = search_point->yy0 * geometry->invstepy + geometry->offsety; + search_point->ip = (int)sx; + search_point->jp = (int)sy; /* test outside of raster */ - if ((ip < 0) || (ip >= n) || (jp < 0) || (jp >= m)) + if ((search_point->ip < 0) || (search_point->ip >= geometry->n) || + (search_point->jp < 0) || (search_point->jp >= geometry->m)) return (3); - if ((ip != iold) || (jp != jold)) { - dx = (double)ip * stepx; - dy = (double)jp * stepy; + if ((search_point->ip != iold) || (search_point->jp != jold)) { + double dx = (double)search_point->ip * geometry->stepx; + double dy = (double)search_point->jp * geometry->stepy; - length = distance( - xg0, dx, yg0, - dy); /* dist from orig. grid point to the current grid point */ - succes2 = test_low_res(); + horizon->length = + distance(origin_point->xg0, dx, origin_point->yg0, dy, + origin_point->coslatsq); /* dist from orig. grid point + to the current grid point */ + int succes2 = test_low_res(geometry, origin_point, origin_angle, + search_point, horizon); if (succes2 == 1) { - zp = z[jp][ip]; + search_point->zp = z[search_point->jp][search_point->ip]; return (1); } } @@ -927,61 +858,68 @@ int new_point(void) return -1; } -int test_low_res(void) +int test_low_res(const Geometry *geometry, const OriginPoint *origin_point, + const OriginAngle *origin_angle, SearchPoint *search_point, + const HorizonProperties *horizon) { - int iold100, jold100; - double sx, sy; - int delx, dely, mindel; - double zp100, z2, curvature_diff; - - iold100 = ip100; - jold100 = jp100; - ip100 = floor(ip / 100.); - jp100 = floor(jp / 100.); + int iold100 = search_point->ip100; + int jold100 = search_point->jp100; + search_point->ip100 = floor(search_point->ip / 100.); + search_point->jp100 = floor(search_point->jp / 100.); /*test the new position with low resolution */ - if ((ip100 != iold100) || (jp100 != jold100)) { - G_debug(2, "ip:%d jp:%d iold100:%d jold100:%d\n", ip, jp, iold100, - jold100); + if ((search_point->ip100 != iold100) || (search_point->jp100 != jold100)) { + G_debug(2, "ip:%d jp:%d iold100:%d jold100:%d\n", search_point->ip, + search_point->jp, iold100, jold100); /* replace with approximate version curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS)); */ - curvature_diff = 0.5 * length * length * invEarth; - z2 = z_orig + curvature_diff + length * tanh0; - zp100 = z100[jp100][ip100]; - G_debug(2, "ip:%d jp:%d z2:%lf zp100:%lf \n", ip, jp, z2, zp100); + double curvature_diff = + 0.5 * horizon->length * horizon->length * invEarth; + double z2 = origin_point->z_orig + curvature_diff + + horizon->length * horizon->tanh0; + double zp100 = z100[search_point->jp100][search_point->ip100]; + G_debug(2, "ip:%d jp:%d z2:%lf zp100:%lf \n", search_point->ip, + search_point->jp, z2, zp100); if (zp100 <= z2) /*skip to the next lowres cell */ { - delx = 32000; - dely = 32000; - if (cosangle > 0.) { - sx = xx0 * invstepx + offsetx; - delx = - floor(fabs((ceil(sx / 100.) - (sx / 100.)) * distcosangle)); + int delx = 32000; + int dely = 32000; + if (origin_angle->cosangle > 0.) { + double sx = + search_point->xx0 * geometry->invstepx + geometry->offsetx; + delx = floor(fabs((ceil(sx / 100.) - (sx / 100.)) * + origin_angle->distcosangle)); } - if (cosangle < 0.) { - sx = xx0 * invstepx + offsetx; - delx = floor( - fabs((floor(sx / 100.) - (sx / 100.)) * distcosangle)); + if (origin_angle->cosangle < 0.) { + double sx = + search_point->xx0 * geometry->invstepx + geometry->offsetx; + delx = floor(fabs((floor(sx / 100.) - (sx / 100.)) * + origin_angle->distcosangle)); } - if (sinangle > 0.) { - sy = yy0 * invstepy + offsety; - dely = - floor(fabs((ceil(sy / 100.) - (sy / 100.)) * distsinangle)); + if (origin_angle->sinangle > 0.) { + double sy = + search_point->yy0 * geometry->invstepy + geometry->offsety; + dely = floor(fabs((ceil(sy / 100.) - (sy / 100.)) * + origin_angle->distsinangle)); } - else if (sinangle < 0.) { - sy = yy0 * invstepy + offsety; - dely = floor( - fabs((floor(sy / 100.) - (sy / 100.)) * distsinangle)); + else if (origin_angle->sinangle < 0.) { + double sy = + search_point->yy0 * geometry->invstepy + geometry->offsety; + dely = floor(fabs((floor(sy / 100.) - (sy / 100.)) * + origin_angle->distsinangle)); } - mindel = min(delx, dely); - G_debug(2, "%d %d %d %lf %lf\n", ip, jp, mindel, xg0, yg0); + int mindel = MIN(delx, dely); + G_debug(2, "%d %d %d %lf %lf\n", search_point->ip, search_point->jp, + mindel, origin_point->xg0, origin_point->yg0); - yy0 = yy0 + (mindel * stepsinangle); - xx0 = xx0 + (mindel * stepcosangle); - G_debug(2, " %lf %lf\n", xx0, yy0); + search_point->yy0 = + search_point->yy0 + (mindel * origin_angle->stepsinangle); + search_point->xx0 = + search_point->xx0 + (mindel * origin_angle->stepcosangle); + G_debug(2, " %lf %lf\n", search_point->xx0, search_point->yy0); return (3); } @@ -995,285 +933,215 @@ int test_low_res(void) } } -double searching(void) +double horizon_height(const Geometry *geometry, const OriginPoint *origin_point, + const OriginAngle *origin_angle) { - double z2; - double curvature_diff; - int succes = 1; + SearchPoint search_point; + HorizonProperties horizon; - if (zp == UNDEFZ) + search_point.ip = 0; + search_point.jp = 0; + search_point.xx0 = origin_point->xg0; + search_point.yy0 = origin_point->yg0; + search_point.zp = origin_point->z_orig; + search_point.ip100 = floor(origin_point->xg0 * geometry->invstepx / 100.); + search_point.jp100 = floor(origin_point->yg0 * geometry->invstepy / 100.); + + horizon.length = 0; + horizon.tanh0 = 0; + + if (search_point.zp == UNDEFZ) return 0; while (1) { - succes = new_point(); + int succes = new_point(geometry, origin_point, origin_angle, + &search_point, &horizon); if (succes != 1) { break; } /* curvature_diff = EARTHRADIUS*(1.-cos(length/EARTHRADIUS)); */ - curvature_diff = 0.5 * length * length * invEarth; + double curvature_diff = + 0.5 * horizon.length * horizon.length * invEarth; - z2 = z_orig + curvature_diff + length * tanh0; + double z2 = origin_point->z_orig + curvature_diff + + horizon.length * horizon.tanh0; - if (z2 < zp) { - tanh0 = (zp - z_orig - curvature_diff) / length; + if (z2 < search_point.zp) { + horizon.tanh0 = + (search_point.zp - origin_point->z_orig - curvature_diff) / + horizon.length; } - if (z2 >= zmax) { + if (z2 >= geometry->zmax) { break; } - if (length >= maxlength) { + if (horizon.length >= origin_point->maxlength) { break; } } - return atan(tanh0); + return atan(horizon.tanh0); } /*////////////////////////////////////////////////////////////////////// */ -void calculate(double xcoord, double ycoord, int buffer_e, int buffer_w, - int buffer_s, int buffer_n) +void calculate_raster_mode(const Settings *settings, const Geometry *geometry, + struct Cell_head *cellhd, + struct Cell_head *new_cellhd, int buffer_e, + int buffer_w, int buffer_s, int buffer_n, + double bufferZone) { - int i, j, l, k = 0; - size_t decimals; - - int xindex, yindex; - double shadow_angle; - double coslat; - - double latitude, longitude; - double xp, yp; - double inputAngle; - double delt_lat, delt_lon; - double delt_east, delt_nor; - double delt_dist; - - char msg_buff[256]; - int hor_row_start = buffer_s; - int hor_row_end = m - buffer_n; + int hor_row_end = geometry->m - buffer_n; int hor_col_start = buffer_w; - int hor_col_end = n - buffer_e; - - int hor_numrows = m - (buffer_s + buffer_n); - int hor_numcols = n - (buffer_e + buffer_w); - - int arrayNumInt; - double dfr_rad, angle_deg; + int hor_col_end = geometry->n - buffer_e; - xindex = (int)((xcoord - xmin) / stepx); - yindex = (int)((ycoord - ymin) / stepy); + int hor_numrows = geometry->m - (buffer_s + buffer_n); + int hor_numcols = geometry->n - (buffer_e + buffer_w); if ((G_projection() == PROJECTION_LL)) { - ll_correction = TRUE; + ll_correction = true; } - if (isMode() == SINGLE_POINT) { - /* Calculate the horizon for one single point */ + /****************************************************************/ + /* The loop over raster points starts here! */ - /* - xg0 = xx0 = (double)xcoord * stepx; - yg0 = yy0 = (double)ycoord * stepy; - xg0 = xx0 = xcoord -0.5*stepx -xmin; - yg0 = yy0 = ycoord -0.5*stepy-ymin; - xg0 = xx0 = xindex*stepx -0.5*stepx; - yg0 = yy0 = yindex*stepy -0.5*stepy; - */ - xg0 = xx0 = xindex * stepx; - yg0 = yy0 = yindex * stepy; + /****************************************************************/ - if (ll_correction) { - coslat = cos(deg2rad * (ymin + yy0)); - coslatsq = coslat * coslat; + if (settings->horizon_basename != NULL) { + horizon_raster = (float **)G_malloc(sizeof(float *) * (hor_numrows)); + for (int l = 0; l < hor_numrows; l++) { + horizon_raster[l] = + (float *)G_malloc(sizeof(float) * (hor_numcols)); } - z_orig = zp = z[yindex][xindex]; - G_debug(1, "yindex: %d, xindex %d, z_orig %.2f", yindex, xindex, - z_orig); - - calculate_shadow(); - fclose(fp); + for (int j = 0; j < hor_numrows; j++) { + for (int i = 0; i < hor_numcols; i++) + horizon_raster[j][i] = 0.; + } + } + double dfr_rad; + int arrayNumInt; + /* definition of horizon angle in loop */ + char *shad_filename = NULL; + if (settings->step == 0.0) { + dfr_rad = 0; + arrayNumInt = 1; + sprintf(shad_filename, "%s", settings->horizon_basename); } else { + dfr_rad = settings->step * deg2rad; + arrayNumInt = 0; + for (double tmp = 0; tmp < settings->end - settings->start; + tmp += fabs(settings->step)) + ++arrayNumInt; + } - /****************************************************************/ - /* The loop over raster points starts here! */ - - /****************************************************************/ - - if (horizon != NULL) { - horizon_raster = - (float **)G_malloc(sizeof(float *) * (hor_numrows)); - for (l = 0; l < hor_numrows; l++) { - horizon_raster[l] = - (float *)G_malloc(sizeof(float) * (hor_numcols)); - } - - for (j = 0; j < hor_numrows; j++) { - for (i = 0; i < hor_numcols; i++) - horizon_raster[j][i] = 0.; - } - } - - /* definition of horizon angle in loop */ - if (step == 0.0) { - dfr_rad = 0; - arrayNumInt = 1; - sprintf(shad_filename, "%s", horizon); - } - else { - dfr_rad = step * deg2rad; - arrayNumInt = 0; - for (double tmp = 0; tmp < end - start; tmp += fabs(step)) - ++arrayNumInt; - } - - decimals = G_get_num_decimals(str_step); - - for (k = 0; k < arrayNumInt; k++) { - struct History history; - - angle = (start + single_direction) * deg2rad + (dfr_rad * k); - angle_deg = angle * rad2deg + 0.0001; - - if (step != 0.0) - shad_filename = - G_generate_basename(horizon, angle_deg, 3, decimals); - - /* - com_par(angle); - */ - G_message( - _("Calculating map %01d of %01d (angle %.2f, raster map <%s>)"), - (k + 1), arrayNumInt, angle_deg, shad_filename); - - for (j = hor_row_start; j < hor_row_end; j++) { - G_percent(j - hor_row_start, hor_numrows - 1, 2); - shadow_angle = 15 * deg2rad; - for (i = hor_col_start; i < hor_col_end; i++) { - ip100 = floor(i / 100.); - jp100 = floor(j / 100.); - ip = jp = 0; - xg0 = xx0 = (double)i * stepx; - - xp = xmin + xx0; - yg0 = yy0 = (double)j * stepy; - - yp = ymin + yy0; - length = 0; - if (ll_correction) { - coslat = cos(deg2rad * yp); - coslatsq = coslat * coslat; - } - - longitude = xp; - latitude = yp; - - if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD, - &longitude, &latitude, NULL) < 0) - G_fatal_error(_("Error in %s"), "GPJ_transform()"); + size_t decimals = G_get_num_decimals(settings->str_step); + + for (int k = 0; k < arrayNumInt; k++) { + struct History history; + + double angle = + (settings->start + settings->single_direction) * deg2rad + + (dfr_rad * k); + double angle_deg = angle * rad2deg + 0.0001; + + if (settings->step != 0.0) + shad_filename = G_generate_basename(settings->horizon_basename, + angle_deg, 3, decimals); + G_message( + _("Calculating map %01d of %01d (angle %.2f, raster map <%s>)"), + (k + 1), arrayNumInt, angle_deg, shad_filename); + + for (int j = hor_row_start; j < hor_row_end; j++) { + G_percent(j - hor_row_start, hor_numrows - 1, 2); + for (int i = hor_col_start; i < hor_col_end; i++) { + OriginPoint origin_point; + OriginAngle origin_angle; + origin_point.xg0 = (double)i * geometry->stepx; + + double xp = geometry->xmin + origin_point.xg0; + origin_point.yg0 = (double)j * geometry->stepy; + + double yp = geometry->ymin + origin_point.yg0; + origin_point.coslatsq = 0; + if (ll_correction) { + double coslat = cos(deg2rad * yp); + origin_point.coslatsq = coslat * coslat; + } - latitude *= deg2rad; - longitude *= deg2rad; + double inputAngle = angle + pihalf; + inputAngle = + (inputAngle >= twopi) ? inputAngle - twopi : inputAngle; + com_par(geometry, &origin_angle, inputAngle, xp, yp); - inputAngle = angle + pihalf; - inputAngle = - (inputAngle >= twopi) ? inputAngle - twopi : inputAngle; + origin_point.z_orig = z[j][i]; + origin_point.maxlength = + (geometry->zmax - origin_point.z_orig) / TANMINANGLE; + origin_point.maxlength = + (origin_point.maxlength < settings->fixedMaxLength) + ? origin_point.maxlength + : settings->fixedMaxLength; - delt_lat = - -0.0001 * cos(inputAngle); /* Arbitrary small distance - in latitude */ - delt_lon = 0.0001 * sin(inputAngle) / cos(latitude); + if (origin_point.z_orig != UNDEFZ) { - latitude = (latitude + delt_lat) * rad2deg; - longitude = (longitude + delt_lon) * rad2deg; + G_debug(4, "**************new line %d %d\n", i, j); + double shadow_angle = + horizon_height(geometry, &origin_point, &origin_angle); - if ((G_projection() != PROJECTION_LL)) { - if (GPJ_transform(&iproj, &oproj, &tproj, PJ_INV, - &longitude, &latitude, NULL) < 0) - G_fatal_error(_("Error in %s"), "GPJ_transform()"); + if (settings->degreeOutput) { + shadow_angle *= rad2deg; } - delt_east = longitude - xp; - delt_nor = latitude - yp; - - delt_dist = - sqrt(delt_east * delt_east + delt_nor * delt_nor); - - sinangle = delt_nor / delt_dist; - cosangle = delt_east / delt_dist; - com_par(); - - z_orig = zp = z[j][i]; - maxlength = (zmax - z_orig) / TANMINANGLE; - maxlength = (maxlength < fixedMaxLength) ? maxlength - : fixedMaxLength; - - if (z_orig != UNDEFZ) { - - G_debug(4, "**************new line %d %d\n", i, j); - shadow_angle = horizon_height(); + horizon_raster[j - buffer_s][i - buffer_w] = shadow_angle; - if (degreeOutput) { - shadow_angle *= rad2deg; - } - - /* - if((j==1400)&&(i==1400)) - { - G_debug(1, "coordinates=%f,%f, raster no. %d, - horizon=%f\n", xp, yp, k, shadow_angle); - } - */ - horizon_raster[j - buffer_s][i - buffer_w] = - shadow_angle; - - } /* undefs */ - } + } /* undefs */ } + } - G_debug(1, "OUTGR() starts..."); - OUTGR(cellhd.rows, cellhd.cols); + G_debug(1, "OUTGR() starts..."); + OUTGR(settings, shad_filename, cellhd); - /* empty array */ - for (j = 0; j < hor_numrows; j++) { - for (i = 0; i < hor_numcols; i++) - horizon_raster[j][i] = 0.; - } + /* empty array */ + for (int j = 0; j < hor_numrows; j++) { + for (int i = 0; i < hor_numcols; i++) + horizon_raster[j][i] = 0.; + } - /* return back the buffered region */ - if (bufferZone > 0.) - Rast_set_window(&new_cellhd); + /* return back the buffered region */ + if (bufferZone > 0.) + Rast_set_window(new_cellhd); - /* write metadata */ - Rast_short_history(shad_filename, "raster", &history); + /* write metadata */ + Rast_short_history(shad_filename, "raster", &history); - sprintf(msg_buff, - "Angular height of terrain horizon, map %01d of %01d", - (k + 1), arrayNumInt); - Rast_put_cell_title(shad_filename, msg_buff); + char msg_buff[256]; + sprintf(msg_buff, "Angular height of terrain horizon, map %01d of %01d", + (k + 1), arrayNumInt); + Rast_put_cell_title(shad_filename, msg_buff); - if (degreeOutput) - Rast_write_units(shad_filename, "degrees"); - else - Rast_write_units(shad_filename, "radians"); + if (settings->degreeOutput) + Rast_write_units(shad_filename, "degrees"); + else + Rast_write_units(shad_filename, "radians"); - Rast_command_history(&history); + Rast_command_history(&history); - /* insert a blank line */ - Rast_append_history(&history, ""); + /* insert a blank line */ + Rast_append_history(&history, ""); - Rast_append_format_history( - &history, - "Horizon view from azimuth angle %.2f degrees CCW from East", - angle * rad2deg); + Rast_append_format_history( + &history, + "Horizon view from azimuth angle %.2f degrees CCW from East", + angle * rad2deg); - Rast_write_history(shad_filename, &history); + Rast_write_history(shad_filename, &history); + if (shad_filename) G_free(shad_filename); - } } } diff --git a/raster/r.horizon/testsuite/test_r_horizon.py b/raster/r.horizon/testsuite/test_r_horizon.py index 07a1398b6e9..a085bfb00a0 100644 --- a/raster/r.horizon/testsuite/test_r_horizon.py +++ b/raster/r.horizon/testsuite/test_r_horizon.py @@ -5,7 +5,7 @@ PURPOSE: Test r.horizon -COPYRIGHT: (C) 2015 Anna Petrasova +COPYRIGHT: (C) 2015-2024 Anna Petrasova This program is free software under the GNU General Public License (>=v2). Read the file COPYING that comes with GRASS @@ -15,6 +15,7 @@ from grass.gunittest.case import TestCase from grass.gunittest.main import test from grass.gunittest.gmodules import SimpleModule +from grass.script import raster_what ref1 = """azimuth,horizon_height @@ -109,6 +110,9 @@ def tearDownClass(cls): cls.runModule("g.remove", flags="f", type="raster", name=cls.circle) cls.del_temp_region() + def setUp(self): + self.runModule("g.region", raster="elevation") + def tearDown(self): """Remove horizon map after each test method""" self.runModule("g.remove", flags="f", type="raster", name=self.horizon) @@ -165,7 +169,7 @@ def test_point_mode_multiple_direction_artificial(self): self.assertMultiLineEqual(first=ref4, second=stdout) def test_raster_mode_one_direction(self): - """Test mode with 1 point and one direction""" + """Test mode with one direction and against point mode""" module = SimpleModule( "r.horizon", elevation="elevation", output=self.horizon_output, direction=50 ) @@ -175,12 +179,38 @@ def test_raster_mode_one_direction(self): "max": 0.70678365230560, "stddev": 0.0360724286360789, } + output = "test_horizon_output_from_elevation_050" self.assertRasterFitsUnivar( - raster="test_horizon_output_from_elevation_050", + raster=output, reference=ref, precision=1e6, ) + # test if point mode matches raster mode + coordinates = [ + (634725, 216185), + (633315, 217595), + (633555, 223405), + (639955, 220605), + (637505, 219705), + (641105, 222225), + ] + for coordinate in coordinates: + module = SimpleModule( + "r.horizon", + elevation="elevation", + coordinates=coordinate, + output=self.horizon, + direction=50, + step=0, + ) + self.assertModule(module) + stdout = module.outputs.stdout + first = float(stdout.splitlines()[-1].split(",")[-1]) + what = raster_what(output, coord=coordinate) + second = float(what[0][output]["value"]) + self.assertAlmostEqual(first=first, second=second, delta=0.000001) + def test_raster_mode_multiple_direction(self): self.runModule("g.region", raster="elevation", res=100) module = SimpleModule(