get_area_len from .osm Polygon元素

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

  

posted @ 2021-06-13 21:27  土博姜山山  阅读(38)  评论(0编辑  收藏  举报