From 2c932490c555ac274611ad5ed0e496a78c983efd Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Mon, 8 Apr 2024 09:29:59 -0400 Subject: [PATCH] r.horizon: add distance into output (#3565) --- raster/r.horizon/main.c | 62 ++++++++++++++----- raster/r.horizon/r.horizon.html | 4 ++ raster/r.horizon/testsuite/test_r_horizon.py | 65 ++++++++++++++++++++ 3 files changed, 116 insertions(+), 15 deletions(-) diff --git a/raster/r.horizon/main.c b/raster/r.horizon/main.c index 2e4aaf190f5..24d96ec873f 100644 --- a/raster/r.horizon/main.c +++ b/raster/r.horizon/main.c @@ -104,6 +104,7 @@ typedef struct { } Geometry; typedef struct { + bool horizonDistance; int degreeOutput; int compassOutput; double fixedMaxLength; @@ -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); @@ -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]); @@ -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 = @@ -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 " @@ -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(); @@ -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; @@ -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; } } @@ -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); @@ -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; @@ -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, @@ -1051,7 +1082,7 @@ double horizon_height(const Geometry *geometry, const OriginPoint *origin_point, } } - return atan(horizon.tanh0); + return horizon; } /*////////////////////////////////////////////////////////////////////// */ @@ -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; diff --git a/raster/r.horizon/r.horizon.html b/raster/r.horizon/r.horizon.html index 418ece9571f..1b9972cc702 100644 --- a/raster/r.horizon/r.horizon.html +++ b/raster/r.horizon/r.horizon.html @@ -30,6 +30,10 @@

DESCRIPTION

Using the -c flag, the azimuthal angles will be printed in compass orientation (North=0, clockwise). +

+Activating the -l flag allows to additionally print the distance +to each horizon angle. +

Input parameters:

The elevation parameter is an input elevation raster map. If the buffer options are used (see below), this raster should extend diff --git a/raster/r.horizon/testsuite/test_r_horizon.py b/raster/r.horizon/testsuite/test_r_horizon.py index 2df976e6370..bc4043372d2 100644 --- a/raster/r.horizon/testsuite/test_r_horizon.py +++ b/raster/r.horizon/testsuite/test_r_horizon.py @@ -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" @@ -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(