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
