北京6环边界Geo范围的外边长与区域面积计算

from shapely.geometry import Polygon
import geopandas as gpd
# 粗略的北京六环道路  见 北京6环边界Geo范围数据 - by BBBike https://www.cnblogs.com/jeshy/p/14868463.html
poly = Polygon([(116.072, 39.714),
        (116.075, 39.705),
        (116.078, 39.695),
        (116.099, 39.688),
        (116.122, 39.688),
        (116.167, 39.685),
        (116.202, 39.684),
        (116.275, 39.682),
        (116.37, 39.701),
        (116.415, 39.711),
        (116.459, 39.715),
        (116.516, 39.729),
        (116.543, 39.734),
        (116.604, 39.749),
        (116.627, 39.769),
        (116.648, 39.788),
        (116.665, 39.806),
        (116.673, 39.815),
        (116.685, 39.824),
        (116.691, 39.832),
        (116.697, 39.841),
        (116.707, 39.864),
        (116.712, 39.874),
        (116.715, 39.884),
        (116.718, 39.904),
        (116.718, 39.922),
        (116.714, 39.938),
        (116.715, 39.952),
        (116.716, 39.965),
        (116.709, 39.993),
        (116.7, 40.005),
        (116.69, 40.016),
        (116.661, 40.057),
        (116.656, 40.064),
        (116.652, 40.069),
        (116.643, 40.078),
        (116.64, 40.082),
        (116.637, 40.086),
        (116.633, 40.091),
        (116.63, 40.103),
        (116.629, 40.107),
        (116.627, 40.112),
        (116.624, 40.121),
        (116.623, 40.132),
        (116.617, 40.142),
        (116.613, 40.145),
        (116.608, 40.146),
        (116.601, 40.147),
        (116.584, 40.149),
        (116.542, 40.152),
        (116.501, 40.161),
        (116.468, 40.164),
        (116.399, 40.164),
        (116.344, 40.178),
        (116.279, 40.181),
        (116.173, 40.175),
        (116.15, 40.161),
        (116.138, 40.126),
        (116.129, 40.086),
        (116.11, 40.066),
        (116.08, 40.029),
        (116.071, 39.98),
        (116.099, 39.94),
        (116.122, 39.909),
        (116.122, 39.885),
        (116.11, 39.859),
        (116.096, 39.814),
        (116.091, 39.797),
        (116.084, 39.787),
        (116.077, 39.75),
        (116.073, 39.734)])

crs = {'init': 'epsg:4326'}
gpoly = gpd.GeoSeries(poly, 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): 2487.234688144363
# length(km): 191.11631267814698

#############################################################

from shapely.geometry import Polygon
import xml.etree.cElementTree as et
import matplotlib.pyplot as plt
import geopandas as gpd
import numpy as np

# 精确的北京六环道路 见 beijing_car_seg_6ring road.osm： （北京六环附近及往内的可驾驶道路网络（路网graph为连通图））https://www.cnblogs.com/jeshy/p/14763489.html
# beijing_car_seg_6ring road.osm as the input data
root = et.parse(r"./beijing_car_seg_6ring road.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')
    node_s, node_e = int(element_nd[0].attrib['ref']), int(element_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(E[0], '...', E[-2], E[-1])
print(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): 2270.593649051131
# length(km): 187.67214535661333

北京6环粗略边界Geo范围数据 - by BBBike