from shapely.geometry import Polygon import xml.etree.cElementTree as et import matplotlib.pyplot as plt import geopandas as gpd import numpy as np # xzg_omh.osm 为一个Polygon元素 从 JSOM软件自主圈划的区域生产 root = et.parse(r"./xzg_omh.osm").getroot() nodes = {} for node in root.findall('node'): nodes[int(node.attrib['id'])] = (float(node.attrib['lon']), float(node.attrib['lat'])) edges = [] for way in root.findall('way'): element_nd = way.findall('nd') nd_num = len(element_nd) for ind_nd in range(nd_num-1): node_s, node_e = int(element_nd[ind_nd].attrib['ref']), int(element_nd[ind_nd+1].attrib['ref']) path = (node_s, node_e) edges.append(path) edge = edges[0] E = [] E.append(edge) edges.remove(edge) point_s, point_e = nodes[E[0][0]], nodes[E[0][1]] Point_lists = [] Point_lists.append(list(point_s)) Point_lists.append(list(point_e)) while len(edges) > 0: (node_f_s, node_f_e) = E[-1] for (node_n_s, node_n_e) in edges: if node_f_e == node_n_s: E.append((node_n_s, node_n_e)) point_f_e = nodes[node_n_e] Point_lists.append(list(point_f_e)) # print((node_n_s, node_n_e)) edges.remove((node_n_s, node_n_e)) break # Points.pop() print('edges: ', E[0], '...', E[-2], E[-1]) print('Points: ', Point_lists[0], '...', Point_lists[-2], Point_lists[-1]) road_line_arr = np.array(Point_lists) # 转换成二维坐标表示 six_ring_plg = Polygon(road_line_arr) # Polygon crs = {'init': 'epsg:4326'} gpoly = gpd.GeoSeries(six_ring_plg, crs=crs) print('gpoly.crs:', gpoly.crs) gpoly = gpoly.to_crs(epsg=2436) print('gpoly.to_crs:', gpoly.crs) poly = gpoly[0] print('area(km*km):', poly.area/1.0e6) print('length(km):', poly.length/1.0e3) # area(km*km): 22.593649051131 # length(km): 187.67214535661333