Coverage for mapology/hull.py: 19.35%

1"""Determine hull curves from point sets.""" 

2 

3import copy 

4import functools 

5from typing import no_type_check 

6 

7from mapology import log 

8 

9THE_HULLS = { 

10 'type': 'FeatureCollection', 

11 'name': 'Prefix Region Convex Hulls', 

12 'crs': { 

13 'type': 'name', 

14 'properties': { 

15 'name': 'urn:ogc:def:crs:OGC:1.3:CRS84', 

16 }, 

17 }, 

18 'features': [], 

19} 

20HULL_TEMPLATE = { 

21 'type': 'Feature', 

22 'id': '', 

23 'properties': {'name': ''}, 

24 'geometry': {'type': 'Polygon', 'coordinates': []}, 

25} 

26 

27 

28@no_type_check 

29def convex(coords): 

30 """Executes scan to return points in counter-clockwise order that are on the convex hull (Graham).""" 

31 turn_left, turn_right, turn_none = 1, -1, 0 # noqa 

32 

33 def compare(a, b): 

34 return float(a > b) - float(a < b) 

35 

36 def turn(p, q, r): 

37 return compare((q[0] - p[0]) * (r[1] - p[1]) - (r[0] - p[0]) * (q[1] - p[1]), 0) 

38 

39 def keep_left(hull, r): 

40 while len(hull) > 1 and turn(hull[-2], hull[-1], r) != turn_left: 

41 hull.pop() 

42 if not len(hull) or hull[-1] != r: 

43 hull.append(r) 

44 return hull 

45 

46 points = sorted(coords) 

47 lower_hull = functools.reduce(keep_left, points, []) 

48 upper_hull = functools.reduce(keep_left, reversed(points), []) 

49 return lower_hull.extend(upper_hull[i] for i in range(1, len(upper_hull) - 1)) or lower_hull 

50 

51 

52@no_type_check 

53def extract_prefix_hull(prefix, region_name, coords): 

54 prefix_hull = copy.deepcopy(HULL_TEMPLATE) 

55 prefix_hull['id'] = prefix 

56 prefix_hull['properties']['name'] = region_name 

57 

58 hull_coords = [[lon, lat] for lat, lon in convex(coords)] 

59 

60 lon_min = min(lon for lon, _ in hull_coords) 

61 lon_max = max(lon for lon, _ in hull_coords) 

62 if lon_min < -120 and +120 < lon_max: 

63 log.info('Problematic region found (convex hull crossing the outer sign change of longitudes): %s' % prefix) 

64 lon_sign = tuple(+1 if lon > 0 else -1 for lon, _ in hull_coords) 

65 lon_neg_count = sum(-datum for datum in lon_sign if datum < 0) 

66 lon_pos_count = sum(datum for datum in lon_sign if datum > 0) 

67 if lon_pos_count and lon_neg_count: # Does the polygon cross the longitude sign change? 

68 log.info('- (%d positive / %d negative) longitudes: %s' % (lon_pos_count, lon_neg_count, prefix)) 

69 patch_me_up = lon_neg_count > lon_pos_count 

70 add_me = +360 if patch_me_up else -360 

71 if patch_me_up: 

72 log.info('- patching upwards (adding 360 degrees to longitude where negative): %s' % prefix) 

73 hull_coords = [[lon + add_me, lat] if lon < 0 else [lon, lat] for lon, lat in hull_coords] 

74 else: 

75 log.info('- patching downwards (adding 360 degrees to longitude where positive): %s' % prefix) 

76 hull_coords = [[lon + add_me, lat] if lon > 0 else [lon, lat] for lon, lat in hull_coords] 

77 

78 if prefix == 'ET': 

79 log.info('Patching a North-Eastern ear onto the hull: %s' % prefix) 

80 

81 @no_type_check 

82 def et_earify(pair): 

83 """Monkey patch the ET region to offer an ear to select outside of ED.""" 

84 lon_ne = 12.27 # 9333333333334, 

85 lat_ne = 53.91 # 8166666666664 

86 magic_lon = 13.648875 

87 magic_lat = 54.551359 

88 lon, lat = pair 

89 return [magic_lon, magic_lat] if lon > lon_ne and lat > lat_ne else [lon, lat] 

90 

91 hull_coords = [et_earify(pair) for pair in hull_coords] 

92 

93 # problem regions A1, NZ, NF, PA, UH, 

94 prefix_hull['geometry']['coordinates'].append(hull_coords) 

95 return prefix_hull 

96 

97 

98@no_type_check 

99def update_hull_store(hull_store, prefix_hull): 

100 """DRY.""" 

101 hull_store['features'].append(copy.deepcopy(prefix_hull))