Skip to content

Commit

Permalink
r.horizon: add distance into output (#3565)
Browse files Browse the repository at this point in the history
  • Loading branch information
@petrasovaa
petrasovaa committed Apr 8, 2024
1 parent 597e222 commit 2c93249
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 15 deletions.
62 changes: 47 additions & 15 deletions raster/r.horizon/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ typedef struct {
} Geometry;

typedef struct {
bool horizonDistance;
int degreeOutput;
int compassOutput;
double fixedMaxLength;
Expand All @@ -120,8 +121,9 @@ 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);
HorizonProperties 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,
enum OutputFormat format);
Expand Down Expand Up @@ -163,7 +165,7 @@ int main(int argc, char *argv[])
} parm;

struct {
struct Flag *degreeOutput, *compassOutput;
struct Flag *horizonDistance, *degreeOutput, *compassOutput;
} flag;

G_gisinit(argv[0]);
Expand Down Expand Up @@ -312,6 +314,12 @@ int main(int argc, char *argv[])
_("Name of file for output (use output=- for stdout)");
parm.output->guisection = _("Point mode");

flag.horizonDistance = G_define_flag();
flag.horizonDistance->key = 'l';
flag.horizonDistance->description =
_("Include horizon distance in the output");
flag.horizonDistance->guisection = _("Point mode");

flag.degreeOutput = G_define_flag();
flag.degreeOutput->key = 'd';
flag.degreeOutput->description =
Expand Down Expand Up @@ -357,6 +365,7 @@ int main(int argc, char *argv[])
Settings settings;
settings.degreeOutput = flag.degreeOutput->answer;
settings.compassOutput = flag.compassOutput->answer;
settings.horizonDistance = flag.horizonDistance->answer;

if (G_projection() == PROJECTION_LL)
G_important_message(_("Note: In latitude-longitude coordinate system "
Expand Down Expand Up @@ -796,14 +805,18 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,

origin_point.maxlength = settings->fixedMaxLength;
/* JSON variables and formating */
JSON_Value *root_value, *origin_value, *azimuths_value, *horizons_value;
JSON_Array *coordinates, *azimuths, *horizons;
JSON_Value *root_value, *origin_value, *azimuths_value, *horizons_value,
*distances_value;
JSON_Array *coordinates, *azimuths, *horizons, *distances;
JSON_Object *origin;
json_set_float_serialization_format("%lf");

switch (format) {
case PLAIN:
fprintf(fp, "azimuth,horizon_height\n");
fprintf(fp, "azimuth,horizon_height");
if (settings->horizonDistance)
fprintf(fp, ",horizon_distance");
fprintf(fp, "\n");
break;
case JSON:
root_value = json_value_init_array();
Expand All @@ -816,14 +829,17 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,
azimuths = json_value_get_array(azimuths_value);
horizons_value = json_value_init_array();
horizons = json_value_get_array(horizons_value);
distances_value = json_value_init_array();
distances = json_value_get_array(distances_value);
break;
}
for (int i = 0; i < printCount; i++) {
OriginAngle origin_angle;
com_par(geometry, &origin_angle, angle, xp, yp);

double shadow_angle =
HorizonProperties horizon =
horizon_height(geometry, &origin_point, &origin_angle);
double shadow_angle = atan(horizon.tanh0);

if (settings->degreeOutput) {
shadow_angle *= rad2deg;
Expand All @@ -837,22 +853,32 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,
tmpangle = tmpangle - 360.;
switch (format) {
case PLAIN:
fprintf(fp, "%lf,%lf\n", tmpangle, shadow_angle);
fprintf(fp, "%lf,%lf", tmpangle, shadow_angle);
if (settings->horizonDistance)
fprintf(fp, ",%lf", horizon.length);
fprintf(fp, "\n");
break;
case JSON:
json_array_append_number(azimuths, tmpangle);
json_array_append_number(horizons, shadow_angle);
if (settings->horizonDistance)
json_array_append_number(distances, horizon.length);
break;
}
}
else {
switch (format) {
case PLAIN:
fprintf(fp, "%lf,%lf\n", printangle, shadow_angle);
fprintf(fp, "%lf,%lf", printangle, shadow_angle);
if (settings->horizonDistance)
fprintf(fp, ",%lf", horizon.length);
fprintf(fp, "\n");
break;
case JSON:
json_array_append_number(azimuths, printangle);
json_array_append_number(horizons, shadow_angle);
if (settings->horizonDistance)
json_array_append_number(distances, horizon.length);
break;
}
}
Expand All @@ -874,6 +900,8 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,
if (format == JSON) {
json_object_set_value(origin, "azimuth", azimuths_value);
json_object_set_value(origin, "horizon_height", horizons_value);
if (settings->horizonDistance)
json_object_set_value(origin, "horizon_distance", distances_value);
json_array_append_value(coordinates, origin_value);
char *json_string = json_serialize_to_string_pretty(root_value);
fprintf(fp, "%s\n", json_string);
Expand Down Expand Up @@ -1001,8 +1029,9 @@ int test_low_res(const Geometry *geometry, const OriginPoint *origin_point,
}
}

double horizon_height(const Geometry *geometry, const OriginPoint *origin_point,
const OriginAngle *origin_angle)
HorizonProperties horizon_height(const Geometry *geometry,
const OriginPoint *origin_point,
const OriginAngle *origin_angle)
{
SearchPoint search_point;
HorizonProperties horizon;
Expand All @@ -1018,8 +1047,10 @@ double horizon_height(const Geometry *geometry, const OriginPoint *origin_point,
horizon.length = 0;
horizon.tanh0 = 0;

if (search_point.zp == UNDEFZ)
return 0;
if (search_point.zp == UNDEFZ) {
HorizonProperties h = {0, 0};
return h;
}

while (1) {
int succes = new_point(geometry, origin_point, origin_angle,
Expand Down Expand Up @@ -1051,7 +1082,7 @@ double horizon_height(const Geometry *geometry, const OriginPoint *origin_point,
}
}

return atan(horizon.tanh0);
return horizon;
}

/*////////////////////////////////////////////////////////////////////// */
Expand Down Expand Up @@ -1159,8 +1190,9 @@ void calculate_raster_mode(const Settings *settings, const Geometry *geometry,
if (origin_point.z_orig != UNDEFZ) {

G_debug(4, "**************new line %d %d\n", i, j);
double shadow_angle =
HorizonProperties horizon =
horizon_height(geometry, &origin_point, &origin_angle);
double shadow_angle = atan(horizon.tanh0);

if (settings->degreeOutput) {
shadow_angle *= rad2deg;
Expand Down
4 changes: 4 additions & 0 deletions raster/r.horizon/r.horizon.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,10 @@ <h2>DESCRIPTION</h2>
Using the <b>-c</b> flag, the azimuthal angles will be printed in compass
orientation (North=0, clockwise).

<p>
Activating the <b>-l</b> flag allows to additionally print the distance
to each horizon angle.

<h3>Input parameters:</h3>
<p>The <i>elevation</i> parameter is an input elevation raster map. If
the buffer options are used (see below), this raster should extend
Expand Down
65 changes: 65 additions & 0 deletions raster/r.horizon/testsuite/test_r_horizon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,27 @@
340.000000,0.196863
"""

ref5 = """azimuth,horizon_height,horizon_distance
0.000000,0.197017,5010.039920
20.000000,0.196832,5017.668781
40.000000,0.196875,5017.818251
60.000000,0.196689,5017.220346
80.000000,0.196847,5014.299552
100.000000,0.196645,5019.531851
120.000000,0.196969,5014.957627
140.000000,0.196778,5020.328674
160.000000,0.196863,5013.431958
180.000000,0.197017,5010.039920
200.000000,0.196832,5014.229751
220.000000,0.196875,5011.387034
240.000000,0.196689,5017.220346
260.000000,0.196847,5014.299552
280.000000,0.196645,5019.531851
300.000000,0.196969,5014.957627
320.000000,0.196778,5020.328674
340.000000,0.196863,5013.431958
"""


class TestHorizon(TestCase):
circle = "circle"
Expand Down Expand Up @@ -196,6 +217,50 @@ def test_point_mode_multiple_direction_artificial(self):
stdout = module.outputs.stdout
self.assertMultiLineEqual(first=ref4, second=stdout)

def test_point_mode_multiple_direction_artificial_distance(self):
"""With 1 point, more directions on artificial surface, distance in output"""
module = SimpleModule(
"r.horizon",
elevation=self.circle,
coordinates=(637505, 221755),
output=self.horizon,
direction=0,
step=20,
flags="l",
)
self.assertModule(module)
stdout = module.outputs.stdout
self.assertMultiLineEqual(first=ref5, second=stdout)

module = SimpleModule(
"r.horizon",
elevation=self.circle,
coordinates=(637505, 221755),
output=self.horizon,
direction=0,
step=20,
flags="l",
format="json",
)
self.assertModule(module)
stdout = json.loads(module.outputs.stdout)
azimuths = []
horizons = []
distances = []
reference = {}
for line in ref5.splitlines()[1:]:
azimuth, horizon, distance = line.split(",")
azimuths.append(float(azimuth))
horizons.append(float(horizon))
distances.append(float(distance))
reference["x"] = 637505.0
reference["y"] = 221755.0
reference["azimuth"] = azimuths
reference["horizon_height"] = horizons
reference["horizon_distance"] = distances

self.assertListEqual([reference], stdout)

def test_raster_mode_one_direction(self):
"""Test mode with one direction and against point mode"""
module = SimpleModule(
Expand Down

0 comments on commit 2c93249

Please sign in to comment.