Coverage for brodata / gml.py: 76%

1from shapely.geometry import ( 

2 Point, 

3 LineString, 

4 Polygon, 

5 MultiPoint, 

6 MultiLineString, 

7 MultiPolygon, 

8 GeometryCollection, 

9) 

10 

11 

12# GML 3.2 namespace 

13ns = {"gml": "http://www.opengis.net/gml/3.2"} 

14 

15 

16def parse_poslist(poslist_text, dim=2): 

17 """Convert GML posList string to list of coordinate tuples (2D or 3D).""" 

18 numbers = list(map(float, poslist_text.strip().split())) 

19 if len(numbers) % dim != 0: 

20 raise ValueError(f"Number of coordinates not divisible by dimension {dim}") 

21 return [tuple(numbers[i : i + dim]) for i in range(0, len(numbers), dim)] 

22 

23 

24def polygon_from_gml(polygon_node, dim=2): 

25 """Convert GML 3.2 Polygon or PolygonPatch node to Shapely Polygon.""" 

26 exterior_text = polygon_node.findtext( 

27 ".//gml:exterior/gml:LinearRing/gml:posList", namespaces=ns 

28 ) 

29 if not exterior_text: 

30 raise ValueError("Polygon has no exterior") 

31 exterior = parse_poslist(exterior_text, dim=dim) 

32 

33 interiors = [] 

34 for inner in polygon_node.findall( 

35 ".//gml:interior/gml:LinearRing/gml:posList", namespaces=ns 

36 ): 

37 interiors.append(parse_poslist(inner.text, dim=dim)) 

38 

39 return Polygon(exterior, interiors) 

40 

41 

42def multisurface_from_gml(ms_node, dim=2): 

43 """Convert GML 3.2 MultiSurface node to Shapely MultiPolygon (supports Surface/PolygonPatch).""" 

44 polygons = [] 

45 for member in ms_node.findall(".//gml:surfaceMember", namespaces=ns): 

46 surface_node = member.find(".//gml:Surface", namespaces=ns) 

47 if surface_node is not None: 

48 for patch in surface_node.findall(".//gml:PolygonPatch", namespaces=ns): 

49 polygons.append(polygon_from_gml(patch, dim=dim)) 

50 else: 

51 poly_node = member.find(".//gml:Polygon", namespaces=ns) 

52 if poly_node is not None: 

53 polygons.append(polygon_from_gml(poly_node, dim=dim)) 

54 return MultiPolygon(polygons) 

55 

56 

57def multigeometry_from_gml(mg_node, dim=2): 

58 """Convert GML 3.2 MultiGeometry to Shapely GeometryCollection.""" 

59 geometries = [] 

60 for member in mg_node.findall(".//gml:geometryMember", namespaces=ns): 

61 for child in list(member): 

62 geometries.append(parse_geometry(child, dim=dim)) 

63 return GeometryCollection(geometries) 

64 

65 

66def parse_geometry(node, dim=2): 

67 """Parse any GML 3.2 geometry node to Shapely.""" 

68 tag = node.tag.split("}")[-1] 

69 

70 if tag == "Point": 

71 pos_text = node.findtext(".//gml:pos", namespaces=ns) 

72 coords = tuple(map(float, pos_text.strip().split())) 

73 return Point(coords) 

74 

75 elif tag == "LineString": 

76 pos_list = node.findtext(".//gml:posList", namespaces=ns) 

77 return LineString(parse_poslist(pos_list, dim=dim)) 

78 

79 elif tag in ("Polygon", "PolygonPatch"): 

80 return polygon_from_gml(node, dim=dim) 

81 

82 elif tag == "MultiSurface": 

83 return multisurface_from_gml(node, dim=dim) 

84 

85 elif tag == "MultiPolygon": 

86 polygons = [ 

87 polygon_from_gml(p, dim=dim) 

88 for p in node.findall(".//gml:polygonMember/gml:Polygon", namespaces=ns) 

89 ] 

90 return MultiPolygon(polygons) 

91 

92 elif tag == "MultiLineString": 

93 lines = [ 

94 LineString( 

95 parse_poslist(l.findtext(".//gml:posList", namespaces=ns), dim=dim) 

96 ) 

97 for l in node.findall( 

98 ".//gml:lineStringMember/gml:LineString", namespaces=ns 

99 ) 

100 ] 

101 return MultiLineString(lines) 

102 

103 elif tag == "MultiPoint": 

104 points = [ 

105 Point(tuple(map(float, p.findtext(".//gml:pos", namespaces=ns).split()))) 

106 for p in node.findall(".//gml:pointMember/gml:Point", namespaces=ns) 

107 ] 

108 return MultiPoint(points) 

109 

110 elif tag == "MultiGeometry": 

111 return multigeometry_from_gml(node, dim=dim) 

112 

113 else: 

114 raise NotImplementedError(f"GML type {tag} not supported")