Skip to content

Commit

Permalink
r.horizon: add support for multiple points (#3543)
Browse files Browse the repository at this point in the history
  • Loading branch information
@petrasovaa
petrasovaa committed Apr 9, 2024
1 parent 8e1214c commit afce9ea
Show file tree
Hide file tree
Showing 3 changed files with 113 additions and 40 deletions.
96 changes: 63 additions & 33 deletions raster/r.horizon/main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ HorizonProperties horizon_height(const Geometry *geometry,
const OriginAngle *origin_angle);
void calculate_point_mode(const Settings *settings, const Geometry *geometry,
double xcoord, double ycoord, FILE *fp,
enum OutputFormat format);
enum OutputFormat format, JSON_Object *json_origin);
int new_point(const Geometry *geometry, const OriginPoint *origin_point,
const OriginAngle *origin_angle, SearchPoint *search_point,
HorizonProperties *horizon);
Expand Down Expand Up @@ -155,6 +155,7 @@ double distance(double x1, double x2, double y1, double y2, double coslatsq)
int main(int argc, char *argv[])
{
double xcoord, ycoord;
double *xcoords, *ycoords;
enum OutputFormat format;

struct GModule *module;
Expand Down Expand Up @@ -284,7 +285,8 @@ int main(int argc, char *argv[])

parm.coord = G_define_standard_option(G_OPT_M_COORDS);
parm.coord->description =
_("Coordinate for which you want to calculate the horizon");
_("Coordinate(s) for which you want to calculate the horizon");
parm.coord->multiple = YES;
parm.coord->guisection = _("Point mode");

parm.dist = G_define_option();
Expand Down Expand Up @@ -374,6 +376,7 @@ int main(int argc, char *argv[])
const char *elevin = parm.elevin->answer;
FILE *fp = NULL;
int mode;
int point_count = 0;
if (parm.coord->answer == NULL) {
G_debug(1, "Setting mode: WHOLE_RASTER");
mode = WHOLE_RASTER;
Expand All @@ -385,16 +388,28 @@ int main(int argc, char *argv[])
format = JSON;
else
format = PLAIN;
if (sscanf(parm.coord->answer, "%lf,%lf", &xcoord, &ycoord) != 2) {
G_fatal_error(_(
"Can't read the coordinates from the \"coordinate\" option."));
}

if (xcoord < cellhd.west || xcoord >= cellhd.east ||
ycoord <= cellhd.south || ycoord > cellhd.north) {
G_fatal_error(_("Coordinates are outside of the current region"));
for (int i = 0; parm.coord->answers[i]; i += 2) {
point_count++;
}
xcoords = (double *)G_malloc(point_count * sizeof(double));
ycoords = (double *)G_malloc(point_count * sizeof(double));
for (int i = 0; i < point_count; ++i) {
if (!(G_scan_easting(parm.coord->answers[2 * i], &xcoord,
G_projection()) &&
G_scan_northing(parm.coord->answers[2 * i + 1], &ycoord,
G_projection()))) {
G_fatal_error(_("Can't read the coordinates from the "
"\"coordinate\" option."));
}
if (xcoord < cellhd.west || xcoord >= cellhd.east ||
ycoord <= cellhd.south || ycoord > cellhd.north) {
G_fatal_error(
_("Coordinates are outside of the current region"));
}
xcoords[i] = xcoord;
ycoords[i] = ycoord;
}

/* Transform the coordinates to row/column */

/*
Expand Down Expand Up @@ -584,8 +599,34 @@ int main(int argc, char *argv[])

INPUT(&geometry, elevin);
if (mode == SINGLE_POINT) {
/* Calculate the horizon for one single point */
calculate_point_mode(&settings, &geometry, xcoord, ycoord, fp, format);
JSON_Value *root_value, *origin_value;
JSON_Array *coordinates;
JSON_Object *origin;
if (format == JSON) {
root_value = json_value_init_array();
coordinates = json_value_get_array(root_value);
json_set_float_serialization_format("%lf");
}
for (int i = 0; i < point_count; ++i) {
/* Calculate the horizon for each point */
if (format == JSON) {
origin_value = json_value_init_object();
origin = json_value_get_object(origin_value);
}
calculate_point_mode(&settings, &geometry, xcoords[i], ycoords[i],
fp, format, origin);
if (format == JSON)
json_array_append_value(coordinates, origin_value);
}
if (format == JSON) {
char *json_string = json_serialize_to_string_pretty(root_value);
fprintf(fp, "%s\n", json_string);
json_free_serialized_string(json_string);
json_value_free(root_value);
}
fclose(fp);
G_free(xcoords);
G_free(ycoords);
}
else {
calculate_raster_mode(&settings, &geometry, &cellhd, &new_cellhd,
Expand Down Expand Up @@ -763,7 +804,7 @@ void com_par(const Geometry *geometry, OriginAngle *origin_angle, double angle,

void calculate_point_mode(const Settings *settings, const Geometry *geometry,
double xcoord, double ycoord, FILE *fp,
enum OutputFormat format)
enum OutputFormat format, JSON_Object *json_origin)
{
/*
xg0 = xx0 = (double)xcoord * stepx;
Expand Down Expand Up @@ -805,11 +846,8 @@ 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,
*distances_value;
JSON_Array *coordinates, *azimuths, *horizons, *distances;
JSON_Object *origin;
json_set_float_serialization_format("%lf");
JSON_Value *azimuths_value, *horizons_value, *distances_value;
JSON_Array *azimuths, *horizons, *distances;

switch (format) {
case PLAIN:
Expand All @@ -819,12 +857,9 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,
fprintf(fp, "\n");
break;
case JSON:
root_value = json_value_init_array();
coordinates = json_value_get_array(root_value);
origin_value = json_value_init_object();
origin = json_value_get_object(origin_value);
json_object_set_number(origin, "x", xcoord);
json_object_set_number(origin, "y", ycoord);

json_object_set_number(json_origin, "x", xcoord);
json_object_set_number(json_origin, "y", ycoord);
azimuths_value = json_value_init_array();
azimuths = json_value_get_array(azimuths_value);
horizons_value = json_value_init_array();
Expand Down Expand Up @@ -898,17 +933,12 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,
} /* end of for loop over angles */

if (format == JSON) {
json_object_set_value(origin, "azimuth", azimuths_value);
json_object_set_value(origin, "horizon_height", horizons_value);
json_object_set_value(json_origin, "azimuth", azimuths_value);
json_object_set_value(json_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);
json_free_serialized_string(json_string);
json_value_free(root_value);
json_object_set_value(json_origin, "horizon_distance",
distances_value);
}
fclose(fp);
}

/*////////////////////////////////////////////////////////////////////// */
Expand Down
16 changes: 9 additions & 7 deletions raster/r.horizon/r.horizon.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@ <h2>DESCRIPTION</h2>
outline in one of two modes:

<ul>
<li> single point: as a series of horizon
heights in the specified directions from the given point. The results are
<li> point: as a series of horizon
heights in the specified directions from the given point(s). The results are
written to the stdout.
<li> raster: in this case the output is
one or more raster maps, with each point in a raster giving the horizon
Expand Down Expand Up @@ -82,9 +82,10 @@ <h3>Input parameters:</h3>
contribute significantly to the horizon outline. Note that a viewshed
can be calculated with <em>r.viewshed</em>.

<p>The <i>coordinate</i> parameter takes a pair of easting-northing values
<p>The <i>coordinate</i> parameter takes one or
multiple pairs of easting-northing values
in the current coordinate system and calculates the values of angular
height of the horizon around this point. To achieve the
height of the horizon around each point. To achieve the
consistency of the results, the point coordinate is
aligned to the midpoint of the closest elevation raster cell.

Expand All @@ -106,10 +107,11 @@ <h3>Input parameters:</h3>
horizon raster maps. The raster name of each horizon direction
raster will be constructed as <i>basename_</i>ANGLE, where ANGLE
is the angle in degrees with the direction. If you use <b>r.horizon</b>
in the single point mode this option will be ignored.
in the point mode this option will be ignored.

<p>The <i>file</i> parameter allows saving the resulting horizon
angles in a comma separated ASCII file (single point mode only). If
angles in a comma separated ASCII file with option <b>format=plain</b>
or in a JSON file with option <b>format=json</b> (point mode only). If
you use <b>r.horizon</b> in the raster map mode this option will be ignored.

<p>At the moment the elevation and maximum distance must be measured in meters,
Expand Down Expand Up @@ -147,7 +149,7 @@ <h2>EXAMPLES</h2>

The examples are intended for the North Carolina sample dataset.

<h3>Single point mode</h3>
<h3>Point mode</h3>

<b>Example 1</b>: determine horizon angle in 225 degree direction (output
of horizon angles CCW from East):
Expand Down
41 changes: 41 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 @@ -176,6 +176,20 @@ def test_point_mode_multiple_direction(self):
stdout = module.outputs.stdout
self.assertMultiLineEqual(first=ref2, second=stdout)

def test_point_mode_multiple_points_and_directions(self):
"""Test mode with 2 identical points and multiple directions"""
module = SimpleModule(
"r.horizon",
elevation="elevation",
coordinates=(634720, 216180, 634720, 216180),
output=self.horizon,
direction=180,
step=20,
)
self.assertModule(module)
stdout = module.outputs.stdout
self.assertMultiLineEqual(first=ref2 + ref2, second=stdout)

def test_point_mode_multiple_direction_json(self):
"""Test mode with 1 point and multiple directions with JSON"""
module = SimpleModule(
Expand Down Expand Up @@ -203,6 +217,33 @@ def test_point_mode_multiple_direction_json(self):

self.assertListEqual([reference], stdout)

def test_point_mode_multiple_points_and_directions_json(self):
"""Test mode with 2 identical points and multiple directions with JSON"""
module = SimpleModule(
"r.horizon",
elevation="elevation",
coordinates=(634720, 216180, 634720, 216180),
output=self.horizon,
direction=180,
step=20,
format="json",
)
self.assertModule(module)
stdout = json.loads(module.outputs.stdout)
azimuths = []
horizons = []
reference = {}
for line in ref2.splitlines()[1:]:
azimuth, horizon = line.split(",")
azimuths.append(float(azimuth))
horizons.append(float(horizon))
reference["x"] = 634720.0
reference["y"] = 216180.0
reference["azimuth"] = azimuths
reference["horizon_height"] = horizons

self.assertListEqual([reference, reference], stdout)

def test_point_mode_multiple_direction_artificial(self):
"""Test mode with 1 point and multiple directions with artificial surface"""
module = SimpleModule(
Expand Down

0 comments on commit afce9ea

Please sign in to comment.