Skip to content

Commit

Permalink
small fixes:
Browse files Browse the repository at this point in the history
- fix rounding coordinates, use meter precision for height rounding
- code improvements for improved readability
  • Loading branch information
arbakker committed Oct 5, 2023
1 parent bb2ca30 commit 60f9a5f
Show file tree
Hide file tree
Showing 5 changed files with 47 additions and 42 deletions.
76 changes: 43 additions & 33 deletions src/geodense/lib.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,20 @@
]

DEFAULT_MAX_SEGMENT_LENGTH = 200
DEFAULT_PRECISION_GEOGRAPHIC = 9
DEFAULT_PRECISION_PROJECTED = 4
DEFAULT_PRECISION_DISTANCE = 4
DEFAULT_PRECISION_DEGREES = 9
DEFAULT_PRECISION_METERS = 4

SUPPORTED_FILE_FORMATS = {
"ESRI Shapefile": [".shp"],
"FlatGeobuf": [".fgb"],
"GeoJSON": [".geojson", ".json"],
"GML": [".gml"],
"GPKG": [".gpkg"],
}


SINGLE_LAYER_EXTENSIONS = [".json", ".geojson", ".shp", ".fgb"]

ERROR_MESSAGE_UNSUPPORTED_FILE_FORMAT = "Argument {arg_name} {file_path} is of an unsupported fileformat, see list-formats for list of supported file formats"
THREE_DIMENSIONAL = 3
point_type = tuple[float, ...]
Expand Down Expand Up @@ -169,7 +173,7 @@ def get_cmd_result_message(
)
for i, item in enumerate(report):
ft_index, coordinates_indices = item[0][:1], item[0][1:]
distance = round(item[1], DEFAULT_PRECISION_DISTANCE)
distance = round(item[1], DEFAULT_PRECISION_METERS)
ft_report = f" - features{ft_index}.geometry.segments\
[{', '.join([str(x) for x in coordinates_indices])}], distance: {distance}"
if len(report) - 1 != i:
Expand Down Expand Up @@ -293,7 +297,6 @@ def densify_file( # noqa: PLR0913

max_segment_length = abs(max_segment_length or DEFAULT_MAX_SEGMENT_LENGTH)
layer = _get_valid_layer_name(input_file_path, layer)
single_layer_file_ext = [".json", ".geojson"]

with fiona.open(input_file_path, layer=layer) as src:
profile = src.profile
Expand All @@ -316,18 +319,18 @@ def densify_file( # noqa: PLR0913
)

prec = (
DEFAULT_PRECISION_GEOGRAPHIC
DEFAULT_PRECISION_DEGREES
if _crs_is_geographic(crs)
else DEFAULT_PRECISION_PROJECTED
else DEFAULT_PRECISION_METERS
)
# COORDINATE_PRECISION is only a lco (layer creation option) for OGR GeoJSON driver
fun_kwargs = {"COORDINATE_PRECISION": prec} if is_geojson_driver else {}
is_geojson_driver = profile["driver"] == "GeoJSON"
fun_kwargs: dict = {"COORDINATE_PRECISION": prec} if is_geojson_driver else {}

if output_file_ext not in SINGLE_LAYER_EXTENSIONS:
fun_kwargs["layer"] = layer

with fiona.open(
output_file_path, "w", **profile, layer=layer, **fun_kwargs
) if output_file_ext not in single_layer_file_ext else fiona.open(
output_file_path, "w", **profile, **fun_kwargs
) as dst:
with fiona.open(output_file_path, "w", **profile, **fun_kwargs) as dst:
transformer = _get_transformer(crs, TRANSFORM_CRS)
for i, ft in enumerate(src):
try:
Expand Down Expand Up @@ -467,7 +470,10 @@ def interpolate_geodesic(
tuple(
(
*back_transformer.transform(lon, lat),
(height_a + ((i + 1) * delta_height_per_point)),
round(
(height_a + ((i + 1) * delta_height_per_point)),
DEFAULT_PRECISION_METERS,
),
)
)
for i, (lon, lat) in enumerate(zip(r.lons, r.lats))
Expand Down Expand Up @@ -570,7 +576,7 @@ def _add_vertices_to_line_segment(
max_segment_length: float,
densify_in_projection: bool = False,
) -> int:
"""Adds vertices to line segment in place, and returns number of vertices added.
"""Adds vertices to linestring in place, and returns number of vertices added to linestring.
Args:
ft_linesegment (_type_): line segment to add vertices
Expand All @@ -588,26 +594,31 @@ def _add_vertices_to_line_segment(

prec = _get_coord_precision(transformer)
if not densify_in_projection:
p = _round_line_segment(
interpolate_geodesic(a, b, max_segment_length, transformer), prec
p = list(
[
_round_coordinates(x, prec)
for x in interpolate_geodesic(a, b, max_segment_length, transformer)
]
)
else:
p = _round_line_segment(interpolate_src_proj(a, b, max_segment_length), prec)

result = linestring

result[coord_index] = tuple([round(x, prec) for x in result[coord_index]])
result[coord_index + 1] = tuple([round(x, prec) for x in result[coord_index + 1]])

result[coord_index + 1 : coord_index + 1] = p
p = list(
[
_round_coordinates(x, prec)
for x in interpolate_src_proj(a, b, max_segment_length)
]
)

linestring[coord_index] = _round_coordinates(linestring[coord_index], prec)
linestring[coord_index + 1] = _round_coordinates(linestring[coord_index + 1], prec)
linestring[coord_index + 1 : coord_index + 1] = p
return len(p)


def _round_line_segment(
l_segment: list[point_type], precision: int
) -> list[point_type]:
return list([tuple(round(y, precision) for y in x) for x in l_segment])
def _round_coordinates(coordinates: tuple, position_precision: int) -> tuple:
result = tuple([round(x, position_precision) for x in coordinates[:2]])
if len(coordinates) == THREE_DIMENSIONAL:
result = (*result, round(coordinates[2], DEFAULT_PRECISION_METERS))
return result


def _add_vertices_exceeding_max_segment_length(
Expand Down Expand Up @@ -681,10 +692,9 @@ def _crs_is_geographic(crs_string: str) -> bool:
def _get_coord_precision(transformer: Transformer) -> int:
if transformer.source_crs is None:
raise ValueError("transformer.source_crs is None")
is_geographic = transformer.source_crs.is_geographic
coord_precision = DEFAULT_PRECISION_PROJECTED
if is_geographic:
coord_precision = DEFAULT_PRECISION_GEOGRAPHIC
coord_precision = DEFAULT_PRECISION_METERS
if transformer.source_crs.is_geographic:
coord_precision = DEFAULT_PRECISION_DEGREES
return coord_precision


Expand Down
2 changes: 0 additions & 2 deletions tests/test_cli.py
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,6 @@ def test_cli_densify_cmd(mock_command, tmpdir, test_dir):
):
main()

# TODO: test all CLI options
assert mock_command.call_args.kwargs["input_file"] == in_filepath
assert mock_command.call_args.kwargs["output_file"] == out_filepath
assert mock_command.call_args.kwargs["max_segment_length"] == float(
Expand Down Expand Up @@ -116,7 +115,6 @@ def test_cli_check_density_cmd(mock_command, test_dir):
):
main()

# TODO: test all CLI options
assert mock_command.call_args.kwargs["input_file"] == in_filepath
assert mock_command.call_args.kwargs["max_segment_length"] == float(
max_segment_length
Expand Down
1 change: 0 additions & 1 deletion tests/test_densify.py
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,6 @@ def test_densify_file(test_dir):
assert os.path.exists(out_file)


# TODO: test conversion from geojson to geopackage...
@pytest.mark.parametrize(
("input_file", "output_file", "expectation"),
[
Expand Down
2 changes: 0 additions & 2 deletions tests/test_get_hr_report.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,6 @@
get_cmd_result_message,
)

# TODO: add test to test content of error message


def test_get_empty_hr_report(linestring_feature_multiple_linesegments):
feature = json.loads(linestring_feature_multiple_linesegments)
Expand Down
8 changes: 4 additions & 4 deletions tests/test_interpolate.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
import pytest
from geodense.lib import (
DEFAULT_PRECISION_GEOGRAPHIC,
DEFAULT_PRECISION_PROJECTED,
DEFAULT_PRECISION_DEGREES,
DEFAULT_PRECISION_METERS,
TRANSFORM_CRS,
_add_vertices_exceeding_max_segment_length,
_add_vertices_to_line_segment,
Expand Down Expand Up @@ -45,7 +45,7 @@ def test_interpolate_round_projected():

assert all(
[
str(x)[::-1].find(".") == DEFAULT_PRECISION_PROJECTED
str(x)[::-1].find(".") == DEFAULT_PRECISION_METERS
for p in points_proj
for x in p
]
Expand All @@ -63,7 +63,7 @@ def test_interpolate_round_geographic():

assert all(
[
str(x)[::-1].find(".") == DEFAULT_PRECISION_GEOGRAPHIC
str(x)[::-1].find(".") == DEFAULT_PRECISION_DEGREES
for p in points_geog
for x in p
]
Expand Down

0 comments on commit 60f9a5f

Please sign in to comment.