Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

r.horizon: add distance into output #3565

Merged
merged 3 commits into from
Apr 8, 2024
Merged

r.horizon: add distance into output #3565

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
42831f8
r.horizon: add horizon distance to output
petrasovaa Apr 5, 2024
990150f
add test
petrasovaa Apr 5, 2024
ae1c14f
add documentation
petrasovaa Apr 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 47 additions & 15 deletions raster/r.horizon/main.c
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading