Coverage for src/geodense/lib.py: 94%
345 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-10-11 14:11 +0000
1import json
2import logging
3import math
4import os
5import sys
6import tempfile
7from collections.abc import Callable, Iterable, Sequence
8from enum import Enum
9from functools import partial
10from typing import Literal, TextIO, cast
12from geojson_pydantic import (
13 Feature,
14 GeometryCollection,
15 LineString,
16 MultiLineString,
17 MultiPoint,
18 MultiPolygon,
19 Point,
20 Polygon,
21)
22from geojson_pydantic.geometries import Geometry
23from geojson_pydantic.types import LineStringCoords, Position, Position2D, Position3D
24from pydantic import BaseModel
25from pyproj import CRS
26from shapely import LineString as ShpLineString
27from shapely import Point as ShpPoint
29from geodense.geojson import CrsFeatureCollection
30from geodense.models import DEFAULT_PRECISION_METERS, DenseConfig, GeodenseError
31from geodense.types import (
32 GeojsonCoordinates,
33 GeojsonGeomNoGeomCollection,
34 GeojsonObject,
35 Nested,
36 ReportLineString,
37 T,
38)
40TWO_DIMENSIONAL = 2
41THREE_DIMENSIONAL = 3
42DEFAULT_CRS_2D = "OGC:CRS84"
43DEFAULT_CRS_3D = "OGC:CRS84h"
44SUPPORTED_FILE_FORMATS = {
45 "GeoJSON": [".geojson", ".json"],
46}
48logger = logging.getLogger("geodense")
51class InfValCoordinateError(Exception):
52 pass
55def densify_geojson_object(densify_config: DenseConfig, geojson_obj: GeojsonObject) -> GeojsonObject:
56 validate_geom_type(geojson_obj, "densify")
57 _densify_geometry = partial(densify_geometry, densify_config)
58 return traverse_geojson_geometries(geojson_obj, _densify_geometry)
61def densify_geometry(densify_config: DenseConfig, geometry: GeojsonGeomNoGeomCollection) -> None:
62 _densify_line_segment = partial(densify_line_segment, densify_config)
63 result: GeojsonCoordinates = traverse_linestrings_in_coordinates(geometry.coordinates, _densify_line_segment)
64 geometry.coordinates = result
67def densify_line_segment(
68 densify_config: DenseConfig,
69 coords: GeojsonCoordinates,
70) -> None:
71 linestring = cast(LineStringCoords, coords)
72 added_nodes = 0
73 stop = len(linestring) - 1
74 for i, _ in enumerate(linestring[:stop]):
75 added_nodes += _add_vertices_to_line_segment(linestring, i + added_nodes, densify_config)
78def check_density_geojson_object(densify_config: DenseConfig, geojson_obj: GeojsonObject) -> CrsFeatureCollection:
79 validate_geom_type(geojson_obj, "density-check")
80 _check_density_geometry = partial(check_density_geometry, densify_config)
81 result: Nested[ReportLineString] | None = transform_geojson_geometries(geojson_obj, _check_density_geometry)
82 if result is None:
83 raise ValueError("_check_density_geometry returned None")
85 flat_result: list[ReportLineString] = (
86 list( # filter out None values, these occur when point geometries are part of input
87 filter(lambda x: x is not None, _flatten(result))
88 )
89 )
90 report_fc = _report_line_string_to_geojson(flat_result, ":".join(densify_config.src_crs.to_authority()))
91 return report_fc
94def check_density_geometry(
95 densify_config: DenseConfig,
96 geometry: GeojsonGeomNoGeomCollection,
97) -> Nested[ReportLineString] | None:
98 _check_density_linestring = partial(check_density_linestring, densify_config)
99 return transform_linestrings_in_coordinates(geometry.coordinates, _check_density_linestring)
102def check_density_linestring(
103 densify_config: DenseConfig,
104 linestring: LineStringCoords,
105) -> list[ReportLineString]:
106 result = []
108 for k in range(0, len(linestring) - 1):
109 a: Position = linestring[k]
110 b: Position = linestring[k + 1]
112 a_2d = cast(Position2D, a[0:2])
113 b_2d = cast(Position2D, b[0:2])
115 if densify_config.in_projection:
116 linesegment_dist = _cartesian_distance(a_2d, b_2d)
117 else:
118 if densify_config.src_crs.is_projected: # only convert to basegeographic crs if src_proj is projected
119 transformer = densify_config.transformer
120 if transformer is None:
121 raise GeodenseError("transformer cannot be None when src_crs.is_projected=True")
122 a_t = transformer.transform(*a_2d)
123 b_t = transformer.transform(*b_2d)
124 else: # src_crs is geographic do not transform
125 a_t, b_t = (a_2d, b_2d)
126 g = densify_config.geod
128 _, _, geod_dist = g.inv(*a_t, *b_t, return_back_azimuth=True)
129 if math.isnan(geod_dist):
130 raise GeodenseError(
131 f"unable to calculate geodesic distance, output calculation geodesic distance: {geod_dist}, expected: floating-point number"
132 )
133 linesegment_dist = geod_dist
134 if linesegment_dist > (densify_config.max_segment_length + 0.001):
135 result.append((linesegment_dist, (a, b)))
136 return result
139def _validate_dependent_file_args(
140 input_file_path: str,
141 output_file_path: str | None = None,
142 overwrite: bool = False,
143) -> None:
144 if output_file_path is not None and (input_file_path == output_file_path and input_file_path != "-"):
145 raise GeodenseError(
146 f"input_file and output_file arguments must be different, input_file: {input_file_path}, output_file: {output_file_path}"
147 )
148 if output_file_path is not None and output_file_path != "-":
149 if os.path.exists(output_file_path) and not overwrite:
150 raise GeodenseError(f"output_file {output_file_path} already exists")
151 elif os.path.exists(output_file_path) and overwrite:
152 os.remove(output_file_path)
155def check_density_file( # noqa: PLR0913
156 input_file_path: str,
157 max_segment_length: float,
158 density_check_report_path: str | None = None,
159 src_crs: str | None = None,
160 in_projection: bool = False,
161 overwrite: bool = False,
162) -> tuple[bool, str, int]:
163 if density_check_report_path is None:
164 density_check_report_path = os.path.join(tempfile.mkdtemp(), "check-density-report.json")
166 _validate_dependent_file_args(input_file_path, density_check_report_path, overwrite)
168 with open(input_file_path) if input_file_path != "-" else sys.stdin as src:
169 geojson_obj = textio_to_geojson(src)
170 validate_geom_type(geojson_obj, "check-density")
171 has_3d_coords: Has3D = _has_3d_coordinates(geojson_obj)
172 geojson_src_crs = _get_crs_geojson(geojson_obj, input_file_path, src_crs, has_3d_coords)
173 config = DenseConfig(
174 CRS.from_authority(*geojson_src_crs.split(":")),
175 max_segment_length,
176 in_projection=in_projection,
177 )
178 report_fc = check_density_geojson_object(config, geojson_obj)
180 failed_segment_count = len(report_fc.features)
181 check_status = failed_segment_count == 0
183 if not check_status:
184 with open(density_check_report_path, "w") as f:
185 f.write(report_fc.model_dump_json(indent=4, exclude_none=True))
186 return (check_status, density_check_report_path, len(report_fc.features))
189def _report_line_string_to_geojson(
190 report: list[ReportLineString], src_crs_auth_code: str | None
191) -> CrsFeatureCollection:
192 features: list[Feature] = [
193 Feature(
194 type="Feature",
195 properties={"segment_length": x[0]},
196 geometry=LineString(type="LineString", coordinates=list(x[1])),
197 )
198 for x in report
199 ]
200 result = CrsFeatureCollection(features=features, type="FeatureCollection", name="density-check-report")
201 if src_crs_auth_code is not None:
202 result.set_crs_auth_code(src_crs_auth_code)
203 return result
206def densify_file( # noqa: PLR0913
207 input_file_path: str,
208 output_file_path: str,
209 overwrite: bool = False,
210 max_segment_length: float | None = None,
211 densify_in_projection: bool = False,
212 src_crs: str | None = None,
213) -> None:
214 """_summary_
216 Arguments:
217 input_file_path: assumes file exists otherwise raises FileNotFoundError
218 output_file_path: assumes directory of output_file_path exists
220 Keyword Arguments:
221 layer -- layer name, when no specified and multilayer file, first layer will be used (default: {None})
222 max_segment_length -- max segment length to use for densification (default: {None})
223 densify_in_projection -- user src projection for densification (default: {False})
224 src_crs -- override src crs of input file (default: {None})
226 Raises:
227 ValueError: application errors
228 pyproj.exceptions.CRSError: when crs cannot be found by pyproj
230 """
231 _validate_dependent_file_args(input_file_path, output_file_path, overwrite)
232 src: TextIO
233 with open(input_file_path) if input_file_path != "-" else sys.stdin as src:
234 geojson_obj = textio_to_geojson(src)
235 has_3d_coords: Has3D = _has_3d_coordinates(geojson_obj)
236 geojson_src_crs = _get_crs_geojson(geojson_obj, input_file_path, src_crs, has_3d_coords)
237 config = DenseConfig(
238 CRS.from_authority(*geojson_src_crs.split(":")),
239 max_segment_length,
240 densify_in_projection,
241 )
242 densify_geojson_object(config, geojson_obj)
243 if src_crs is not None and isinstance(geojson_obj, CrsFeatureCollection):
244 geojson_obj.set_crs_auth_code(src_crs)
245 with open(output_file_path, "w") if output_file_path != "-" else sys.stdout as out_f:
246 geojson_obj_model: BaseModel = cast(BaseModel, geojson_obj)
248 out_f.write(geojson_obj_model.model_dump_json(indent=1, exclude_none=True))
251def transform_linestrings_in_coordinates(
252 coordinates: GeojsonCoordinates,
253 callback: Callable[[GeojsonCoordinates], T],
254) -> Nested | T | None:
255 # when point geometry return None, since we only want to operate on linestring-like coordinates sequences
256 if isinstance(coordinates, tuple):
257 return None
258 elif not _is_linestring_geom(coordinates):
259 _self = partial(transform_linestrings_in_coordinates, callback=callback)
260 return list(
261 map(
262 _self,
263 coordinates,
264 )
265 )
266 else:
267 return callback(coordinates)
270def traverse_linestrings_in_coordinates(
271 coordinates: GeojsonCoordinates,
272 callback: Callable[[GeojsonCoordinates], None],
273) -> GeojsonCoordinates:
274 """Differs from transform_linestrings_in_coordinates in that it does not transform the type of coordinates, so coordinates in > coordinates out. transform_linestrings_in_coordinates expects a callback that transforms the linestring coordinates to generic type T."""
275 # maybe do perform mutations on copy
276 if isinstance(coordinates, tuple):
277 return coordinates
278 elif not _is_linestring_geom(coordinates):
279 _self = partial(traverse_linestrings_in_coordinates, callback=callback)
280 list(
281 map(
282 _self,
283 coordinates,
284 )
285 )
286 else:
287 callback(coordinates)
288 return coordinates
291def traverse_geojson_geometries(
292 geojson: Feature | CrsFeatureCollection | Geometry,
293 geometry_callback: Callable[[Geometry], None] | None = None,
294 node_callback: Callable | None = None,
295) -> Feature | CrsFeatureCollection | Geometry:
296 _self: Callable[
297 [Feature | CrsFeatureCollection | Geometry],
298 Feature | CrsFeatureCollection | Geometry,
299 ] = partial(
300 traverse_geojson_geometries,
301 geometry_callback=geometry_callback,
302 node_callback=node_callback,
303 )
305 _geojson = geojson.model_copy(deep=True)
307 if isinstance(geojson, Feature):
308 feature = cast(Feature, geojson)
309 _feature = cast(Feature, _geojson)
310 geom = cast(GeojsonGeomNoGeomCollection, feature.geometry)
311 try:
312 _feature.geometry = _self(geom)
313 except InfValCoordinateError as _:
314 _feature.geometry = None
316 elif isinstance(geojson, CrsFeatureCollection):
317 _feature_collection: CrsFeatureCollection = cast(CrsFeatureCollection, _geojson)
318 _feature_collection.features = list(map(_self, _feature_collection.features))
319 elif isinstance(geojson, GeometryCollection):
320 _geometry_collection = cast(GeometryCollection, _geojson)
321 __self = cast(
322 Callable[
323 [Geometry],
324 Geometry,
325 ],
326 _self,
327 ) # cast to more specific type, to fix mypy error
328 _geometry_collection.geometries = list(map(__self, _geometry_collection.geometries))
329 elif isinstance(
330 geojson,
331 Point | MultiPoint | LineString | MultiLineString | Polygon | MultiPolygon,
332 ):
333 if geometry_callback is not None:
334 _geom = cast(GeojsonGeomNoGeomCollection, _geojson)
335 geometry_callback(_geom)
337 if node_callback:
338 node_callback(_geojson)
340 return _geojson
343def transform_geojson_geometries(
344 geojson: (Feature | CrsFeatureCollection | GeojsonGeomNoGeomCollection | GeometryCollection),
345 geometry_callback: Callable[[GeojsonGeomNoGeomCollection], T],
346) -> Nested | T:
347 _self = partial(transform_geojson_geometries, geometry_callback=geometry_callback)
349 if isinstance(geojson, Feature):
350 feature = cast(Feature, geojson)
351 geom = cast(GeojsonGeomNoGeomCollection, feature.geometry)
352 return _self(geom)
354 elif isinstance(geojson, CrsFeatureCollection):
355 feature_collection: CrsFeatureCollection = cast(CrsFeatureCollection, geojson)
356 return list(map(_self, feature_collection.features))
357 elif isinstance(geojson, GeometryCollection):
358 geometry_collection = cast(GeometryCollection, geojson)
359 return list(map(_self, geometry_collection.geometries))
360 else:
361 # isinstance(
362 # geojson,
363 # Point | MultiPoint | LineString | MultiLineString | Polygon | MultiPolygon,
364 # ):
365 geom = cast(GeojsonGeomNoGeomCollection, geojson)
366 return geometry_callback(geom)
369def interpolate_geodesic(a: Position, b: Position, densify_config: DenseConfig) -> LineStringCoords:
370 """geodesic interpolate intermediate points between points a and b, with segment_length < max_segment_length. Only returns intermediate points."""
372 three_dimensional_points = len(a) == THREE_DIMENSIONAL and len(b) == THREE_DIMENSIONAL
373 a_2d = Position2D(longitude=a.longitude, latitude=a.latitude)
374 b_2d = Position2D(longitude=b.longitude, latitude=b.latitude)
376 transformer = densify_config.transformer
378 if densify_config.src_crs.is_projected: # only convert to basegeographic crs if src_proj is projected
379 if transformer is None:
380 raise GeodenseError("transformer cannot be None when src_crs.is_projected=True")
381 # technically the following transform call is a converion and not a transformation, since crs->base-crs will be a conversion in most cases
382 a_t: tuple[float, float] = transformer.transform(*a_2d)
383 b_t: tuple[float, float] = transformer.transform(*b_2d)
384 else: # src_crs is geographic do not transform
385 a_t, b_t = (a_2d, b_2d)
387 g = densify_config.geod
389 az12, _, geod_dist = g.inv(*a_t, *b_t, return_back_azimuth=True)
390 if math.isnan(geod_dist):
391 raise GeodenseError(
392 f"unable to calculate geodesic distance, output calculation geodesic distance: {geod_dist}, expected: floating-point number"
393 )
395 if geod_dist <= densify_config.max_segment_length:
396 return []
397 else:
398 (
399 nr_points,
400 new_max_segment_length,
401 ) = _get_intermediate_nr_points_and_segment_length(geod_dist, densify_config.max_segment_length)
402 r = g.fwd_intermediate(
403 *a_t,
404 az12,
405 npts=nr_points,
406 del_s=new_max_segment_length,
407 return_back_azimuth=True,
408 )
410 def optional_back_transform(lon: float, lat: float) -> Position2D:
411 """technically should be named optional_back_convert, since crs->base crs is (mostly) a conversion and not a transformation"""
412 if densify_config.src_crs.is_projected:
413 if densify_config.back_transformer is None:
414 raise GeodenseError("back_transformer cannot be None when src_crs.is_projected=True")
415 _result = densify_config.back_transformer.transform(lon, lat)
416 return Position2D(longitude=_result[0], latitude=_result[1])
417 return Position2D(longitude=lon, latitude=lat)
419 if three_dimensional_points:
420 # interpolate height for three_dimensional_points
421 a_3d = cast(Position3D, a)
422 b_3d = cast(Position3D, b)
423 height_a = a_3d[2:][0]
424 height_b = b_3d[2:][0]
425 delta_height_b_a = height_b - height_a
426 delta_height_per_point = delta_height_b_a * (new_max_segment_length / geod_dist)
427 return [
428 Position3D(
429 *optional_back_transform(lon, lat),
430 altitude=round(
431 (height_a + ((i + 1) * delta_height_per_point)),
432 DEFAULT_PRECISION_METERS,
433 ),
434 )
435 for i, (lon, lat) in enumerate(zip(r.lons, r.lats, strict=True))
436 ]
437 else:
438 return [optional_back_transform(lon, lat) for lon, lat in zip(r.lons, r.lats, strict=True)]
441def _cartesian_distance(a: Position, b: Position) -> float:
442 return math.sqrt((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2) # pythagoras
445def interpolate_src_proj(a: Position, b: Position, densify_config: DenseConfig) -> LineStringCoords:
446 """Interpolate intermediate points between points a and b, with segment_length < max_segment_length. Only returns intermediate points."""
448 all_three_dimensional = len(a) == THREE_DIMENSIONAL and len(b) == THREE_DIMENSIONAL
449 # when mixed 2D/3D reduce to 2D, shapely cannot do interpolation on points with mixed dimensionality
450 if not all_three_dimensional:
451 a = Position2D(longitude=a.longitude, latitude=a.latitude)
452 b = Position2D(longitude=b.longitude, latitude=b.latitude)
454 dist = _cartesian_distance(a, b)
455 if dist <= densify_config.max_segment_length:
456 return []
457 else:
458 new_points: list[Position] = []
460 (
461 nr_points,
462 new_max_segment_length,
463 ) = _get_intermediate_nr_points_and_segment_length(dist, densify_config.max_segment_length)
465 for i in range(0, nr_points):
466 p_point: ShpPoint = ShpLineString([a, b]).interpolate(new_max_segment_length * (i + 1))
467 p_2d = Position2D(longitude=p_point.coords[0][0], latitude=p_point.coords[0][1])
468 if len(p_point.coords[0]) == THREE_DIMENSIONAL:
469 p: Position = Position3D(*p_2d, altitude=p_point.coords[0][2]) # type: ignore
470 else:
471 p = p_2d
472 new_points.append(p)
473 return [
474 *new_points,
475 ]
478def textio_to_geojson(src: TextIO) -> GeojsonObject:
479 src_json = json.loads(src.read())
480 type_map = {
481 "Feature": Feature,
482 "GeometryCollection": GeometryCollection,
483 "FeatureCollection": CrsFeatureCollection,
484 "Point": Point,
485 "MultiPoint": MultiPoint,
486 "Polygon": Polygon,
487 "MultiPolygon": MultiPolygon,
488 "LineString": LineString,
489 "MultiLineString": MultiLineString,
490 }
491 try:
492 geojson_type = src_json["type"]
493 constructor = type_map[geojson_type]
494 except KeyError as e:
495 message = f'received invalid GeoJSON file, loc: `.type`, value: `{src_json["type"]}`, expected one of: {", ".join(list(type_map.keys()))}'
496 raise GeodenseError(message) from e
497 geojson_obj: GeojsonObject = constructor(**src_json)
498 return geojson_obj
501class Has3D(Enum):
502 all: Literal["all"] = "all"
503 some: Literal["some"] = "some"
504 none: Literal["none"] = "none"
507def _get_crs_geojson(
508 geojson_object: GeojsonObject,
509 input_file_path: str,
510 src_crs: str | None,
511 has_3d_coords: Has3D,
512) -> str:
513 result_crs: str | None = None
515 is_fc = False
516 if isinstance(geojson_object, CrsFeatureCollection):
517 result_crs = geojson_object.get_crs_auth_code()
518 is_fc = True
520 # case to set default CRS if src_crs not specified and not available in GeoJSON
521 if (
522 result_crs is None and src_crs is None
523 ): # set default crs if not in geojson object and not overridden with src_crs
524 default_crs = DEFAULT_CRS_2D
525 if has_3d_coords in [Has3D.some, Has3D.all]:
526 default_crs = DEFAULT_CRS_3D
527 message = f"unable to determine source CRS for file {input_file_path}, assumed CRS is {default_crs}"
528 logger.warning(message)
529 result_crs = default_crs
530 if is_fc:
531 fc: CrsFeatureCollection = cast(CrsFeatureCollection, geojson_object)
532 fc.set_crs_auth_code(result_crs)
534 # is src_crs is set use src_crs
535 elif src_crs is not None:
536 src_crs_crs: CRS = CRS.from_authority(*src_crs.split(":"))
537 if has_3d_coords == Has3D.all and not src_crs_crs.is_vertical:
538 logger.warning("src_crs is 2D while input data contains geometries with 3D coordinates")
539 result_crs = src_crs if src_crs is not None else result_crs # override json_crs with src_crs if defined
541 elif result_crs is None:
542 raise GeodenseError("could not determine CRS from GeoJSON object")
544 return result_crs
547def _flatten(container: Nested) -> Iterable:
548 if isinstance(container, tuple | str):
549 raise ValueError("var container should not be of type tuple or str")
551 for item in container:
552 if isinstance(item, Sequence) and not isinstance(item, tuple) and not isinstance(item, str):
553 yield from _flatten(item)
554 else:
555 yield item
558def _is_linestring_geom(geometry_coordinates: GeojsonCoordinates) -> bool:
559 """Check if coordinates are of linestring geometry type.
561 - Fiona linestring coordinates are of type: list[tuple[float,float,...]])
562 - GeoJSON linestring coordinates are of type: list[list[float]]
564 Args:
565 geometry_coordinates (list): Fiona or GeoJSON coordinates sequence
567 Returns:
568 bool: if geometry_coordinates is linestring return True else False
569 """
571 return (
572 len(geometry_coordinates) > 0
573 and isinstance(geometry_coordinates[0], Sequence)
574 and all(isinstance(x, float | int) for x in geometry_coordinates[0])
575 ) # also test for int just in case...
578def _transform_positions_in_coordinates(
579 coordinates: GeojsonCoordinates,
580 callback: Callable[[GeojsonCoordinates], T],
581) -> Nested | T:
582 if not isinstance(coordinates, tuple):
583 return list(map(lambda x: _transform_positions_in_coordinates(x, callback), coordinates))
584 return callback(coordinates)
587def _get_intermediate_nr_points_and_segment_length(dist: float, max_segment_length: float) -> tuple[int, float]:
588 if dist <= max_segment_length:
589 raise GeodenseError(f"max_segment_length ({max_segment_length}) cannot be bigger or equal than dist ({dist})")
590 remainder = dist % max_segment_length
591 nr_segments = int(dist // max_segment_length)
592 if remainder > 0:
593 nr_segments += 1
594 new_max_segment_length = dist / nr_segments # space segments evenly over delta(a,b)
595 nr_points = nr_segments - 1 # convert nr of segments to nr of intermediate points, should be at least 1
596 return nr_points, new_max_segment_length
599def _add_vertices_to_line_segment(linestring: LineStringCoords, coord_index: int, densify_config: DenseConfig) -> int:
600 """Adds vertices to linestring in place, and returns number of vertices added to linestring.
602 Args:
603 ft_linesegment (_type_): line segment to add vertices
604 coord_index (int): coordinate index of line segment to add vertices for
605 transformer (Transformer): pyproj transformer
606 max_segment_length (float): max segment length, if exceeded vertices will be added
607 densify_in_projection (bool): whether to use source projection to densify (not use great-circle distance)
609 Returns:
610 int: number of added vertices
611 """
613 a = linestring[coord_index]
614 b = linestring[coord_index + 1]
616 prec = densify_config.get_coord_precision()
618 if not densify_config.in_projection:
619 linestring_coords = interpolate_geodesic(a, b, densify_config)
620 else:
621 linestring_coords = interpolate_src_proj(a, b, densify_config)
623 p = list(map(lambda x: _round_coordinates(x, prec), linestring_coords))
625 linestring[coord_index] = _round_coordinates(linestring[coord_index], prec)
626 linestring[coord_index + 1] = _round_coordinates(linestring[coord_index + 1], prec)
627 linestring[coord_index + 1 : coord_index + 1] = p
628 return len(p)
631def _round_coordinates(position: Position, precision: int) -> Position:
632 result: Position = Position2D(
633 longitude=round(position.longitude, precision),
634 latitude=round(position.latitude, precision),
635 )
637 if len(position) == THREE_DIMENSIONAL:
638 position_3d = cast(Position3D, position)
639 result = Position3D(
640 longitude=result.longitude,
641 latitude=result.latitude,
642 altitude=round(position_3d.altitude, DEFAULT_PRECISION_METERS),
643 )
644 return result
647def _get_geometry_type(
648 geometry: GeojsonGeomNoGeomCollection,
649) -> str:
650 return cast(str, geometry.type)
653def _geom_has_3d_coords(
654 geometry: GeojsonGeomNoGeomCollection,
655) -> Nested[bool] | bool:
656 def _position_is_3d(position: GeojsonCoordinates) -> bool:
657 return len(position) == THREE_DIMENSIONAL
659 return _transform_positions_in_coordinates(geometry.coordinates, _position_is_3d)
662def _has_3d_coordinates(geojson_obj: GeojsonObject, silent: bool | None = False) -> Has3D:
663 has_3d_coords: Nested[bool] | bool = transform_geojson_geometries(geojson_obj, _geom_has_3d_coords)
665 # in case only value returned from transform_geojson_geometries
666 if isinstance(has_3d_coords, bool): # noqa: SIM108
667 has_3d_coords_flat = [has_3d_coords]
668 else:
669 has_3d_coords_flat = list(_flatten(has_3d_coords))
671 if not all(has_3d_coords_flat) and any(has_3d_coords_flat): # some 3d
672 if not silent:
673 warning_message = "geometries with mixed 2D and 3D vertices found"
674 logger.warning(warning_message)
675 return Has3D.some
676 elif all(not x for x in has_3d_coords_flat): # none 3d
677 return Has3D.none
678 return Has3D.all
681def validate_geom_type(geojson_obj: GeojsonObject, command: str = "") -> None:
682 geom_types: Nested[str] | str = transform_geojson_geometries(geojson_obj, _get_geometry_type)
683 geom_types = [geom_types] if isinstance(geom_types, str) else geom_types
684 geom_types_flat = list(_flatten(geom_types))
685 if all(g_t in ("Point", "MultiPoint") for g_t in geom_types_flat):
686 # situation: all geoms point -> error
687 if command:
688 error_message = f"cannot run {command} on GeoJSON that only contains (Multi)Point geometries"
689 else:
690 error_message = "GeoJSON contains only (Multi)Point geometries"
691 raise GeodenseError(error_message)
692 elif any(gt in ["Point", "MultiPoint"] for gt in geom_types_flat):
693 # sitation: some geoms point -> warning
694 warning_message = "GeoJSON contains (Multi)Point geometries"
695 if command:
696 warning_message = f"{warning_message}, cannot run {command} on (Multi)Point geometries"
698 logger.warning(warning_message)
699 else:
700 # situation: no geoms point -> ok
701 pass