Coverage for pyguymer3/geo/cleanSrc/clean_CoordinateSequence.py: 66%

44 statements  

« prev     ^ index     » next       coverage.py v7.9.2, created at 2025-07-08 18:47 +0000

1#!/usr/bin/env python3 

2 

3# Define function ... 

4def clean_CoordinateSequence( 

5 coords, 

6 /, 

7 *, 

8 debug = __debug__, 

9 prefix = ".", 

10 tol = 1.0e-10, 

11): 

12 """Clean a CoordinateSequence 

13 

14 This function cleans a CoordinateSequence by removing bad points. 

15 

16 Parameters 

17 ---------- 

18 coords : shapely.coords.CoordinateSequence 

19 the CoordinateSequence 

20 debug : bool, optional 

21 print debug messages 

22 prefix : str, optional 

23 change the name of the output debugging CSVs 

24 tol : float, optional 

25 the Euclidean distance that defines two points as being the same (in 

26 degrees) 

27 

28 Returns 

29 ------- 

30 cleans : shapely.geometry.linestring.LineString 

31 the cleaned CoordinateSequence 

32 

33 Notes 

34 ----- 

35 According to the `Shapely documentation for the function 

36 shapely.geometry.polygon.orient() 

37 <https://shapely.readthedocs.io/en/stable/manual.html#shapely.geometry.polygon.orient>`_ : 

38 

39 "A sign of 1.0 means that the coordinates of the product's exterior ring 

40 will be oriented counter-clockwise." 

41 

42 Copyright 2017 Thomas Guymer [1]_ 

43 

44 References 

45 ---------- 

46 .. [1] PyGuymer3, https://github.com/Guymer/PyGuymer3 

47 """ 

48 

49 # Import special modules ... 

50 try: 

51 import numpy 

52 except: 

53 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None 

54 try: 

55 import shapely 

56 import shapely.geometry 

57 except: 

58 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None 

59 

60 # Import sub-functions ... 

61 from ..check import check 

62 

63 # ************************************************************************** 

64 

65 # Check argument ... 

66 assert isinstance(coords, shapely.coords.CoordinateSequence), "\"coords\" is not a CoordinateSequence" 

67 

68 # Convert the CoordinateSequence to a NumPy array ... 

69 points1 = numpy.array(coords) # [°] 

70 

71 # Create short-hand ... 

72 npoint = int(points1.shape[0]) # [#] 

73 

74 # Don't clean points ... 

75 if npoint == 1: 

76 return shapely.geometry.point.Point(coords) 

77 

78 # Initialize list ... 

79 points2 = [] # [°] 

80 points2.append(points1[0, :]) # [°] 

81 

82 # Initialize flag ... 

83 skip = False 

84 

85 # Loop over points ... 

86 for ipoint in range(1, npoint - 1): 

87 # Skip this point if it needs skipping ... 

88 if skip: 

89 # Re-set flag ... 

90 skip = False 

91 

92 continue 

93 

94 # Find the Euclidean distance between the points either side of this 

95 # point ... 

96 dx = points1[ipoint - 1, 0] - points1[ipoint + 1, 0] # [°] 

97 dy = points1[ipoint - 1, 1] - points1[ipoint + 1, 1] # [°] 

98 dr = numpy.hypot(dx, dy) # [°] 

99 

100 # Skip this point (and the next one) if it is a bad spike ... 

101 if dr < tol: 

102 if debug: 

103 print(f"INFO: Removing spike between ({points1[ipoint - 1, 0]:+.6f}°,{points1[ipoint - 1, 1]:+.6f}°) and ({points1[ipoint + 1, 0]:+.6f}°,{points1[ipoint + 1, 1]:+.6f}°), which are {dr:+.6f}° apart.") 

104 

105 # Set flag ... 

106 skip = True 

107 

108 continue 

109 

110 # Append point to list ... 

111 points2.append(points1[ipoint, :]) # [°] 

112 

113 # Append last point to list (if it shouldn't be skipped) ... 

114 if not skip: 

115 points2.append(points1[-1, :]) # [°] 

116 

117 # Convert list of points into a LineString, or a Point if there ends up 

118 # being only one ... 

119 if len(points2) == 1: 

120 cleans = shapely.geometry.point.Point(points2[0][0], points2[0][1]) 

121 else: 

122 cleans = shapely.geometry.linestring.LineString(points2) 

123 if cleans.length < tol: 

124 if debug: 

125 print(f"INFO: Truncating a tiny-length line at ({points2[0][0]:+.6f}°,{points2[0][1]:+.6f}°).") 

126 cleans = shapely.geometry.point.Point(points2[0][0], points2[0][1]) 

127 if debug: 

128 check(cleans, prefix = prefix) 

129 

130 # Return cleaned CoordinateSequence ... 

131 return cleans