from shapely.geometry import Polygon, mapping
import geopandas as gpd
from matplotlib import pyplot as plt
import pandas as pd
polys1 = gpd.GeoSeries([Polygon([(0,0), (2,0), (2,2), (0,2)]),
Polygon([(2,2), (4,2), (4,4), (2,4)])])
polys2 = gpd.GeoSeries([Polygon([(1,1), (3,1), (3,3), (1,3)]),
Polygon([(3,3), (5,3), (5,5), (3,5)])])
df1 = gpd.GeoDataFrame({'geometry': polys1, 'df1':[1,2]})
df2 = gpd.GeoDataFrame({'geometry': polys2, 'df2':[1,2]})
res_union = gpd.overlay(df1, df2, how='intersection')
# >>> res_union
# >>> df1 df2 geometry
# >>> 0 1 1 POLYGON ((1 2, 2 2, 2 1, 1 1, 1 2))
# >>> 1 2 1 POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))
# >>> 2 2 2 POLYGON ((3 4, 4 4, 4 3, 3 3, 3 4))
x,y = res_union['geometry'][0].exterior.coords.xy
# >>>x,y
# >>>(array('d', [1.0, 2.0, 2.0, 1.0, 1.0]), array('d', [2.0, 2.0, 1.0, 1.0, 2.0]))
coordinates = mapping(res_union['geometry'][0])['coordinates'][0]
# >>> ((1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0), (1.0, 2.0))
ax = df1.plot(color='red')
df2.plot(ax=ax, color='green', alpha=0.5)
plt.show()
res_union.plot(ax=ax,color='yellow',alpha=0.5)
plt.show()
Geopandas 依赖库
从 http://www.lfd.uci.edu/~gohlke/pythonlibs 下载 Fiona , GDAl , pyproj , Shapely,Rtree，descartes
这是两组矩形图，重叠部分可以从颜色上看出来
黄色是后续计算出重叠部分，然后再重新画上颜色
20200327增加
将polygon转成list的小功能
def polygon2list(polys):
'''
shapely空间函数数据点提取
:param polys: shapely.geometry.polygon.Polygon
:return: list
'''
temp_str = polys.to_wkt()
temp_str = temp_str.replace('POLYGON ','').replace('(','').replace(')','')
temp_str = temp_str.split(',')
temp_str = [i.strip().split(' ') for i in temp_str]
temp_str = [[round(float(i[0]),6),round(float(i[1]),6)] for i in temp_str]
return temp_str
额外的画图
# 合并画图
polys = gpd.GeoSeries(source['polygon'].tolist())
b = polys.boundary
c = polys.convex_hull
d = gpd.GeoSeries(polys.unary_union)
e = gpd.GeoSeries(d).boundary
b.plot(alpha=0.5)
plt.savefig(f'pic/all_boundary.jpg')
c.plot(alpha=0.5)
plt.savefig(f'pic/all_convex_hull.jpg')
d.plot(alpha=0.5)
plt.savefig(f'pic/all_unary_union.jpg')
e.plot(alpha=0.5)
plt.savefig(f'pic/all_unary_union_boundary.jpg')
boundary
convex_hull
unary_union
unary_union_boundary
# 20200408补充分块画图
rssult_df['border_polygon'] = rssult_df['border'].apply(lambda x:Polygon(x))
polys = gpd.GeoSeries(rssult_df['border_polygon'].tolist())
b = polys.boundary
c = polys.convex_hull
d = gpd.GeoSeries(polys.unary_union)
e = gpd.GeoSeries(d).boundary
b.plot(alpha=0.5,figsize=(30,20))
plt.savefig(f'pic/all_boundary_1.jpg')
20200706增加
使用得是shapely来做计算
from shapely.geometry import Point,MultiPoint,Polygon
def createArea(lon_lat:list,area_tyoe = 'rectangle', radius = 0.009):
'''
创建区域(标准方形还是圆形)
:param lon_lat: 经纬度
:param area_tyoe: 创建类型:rectangle 方形 roundness 圆形
:param radius: 方形为边长1半，圆形为半径
:return: 图形边界数据
'''
if area_tyoe == 'rectangle':
return Polygon([[lon_lat[0] + radius, lon_lat[1] + radius],
[lon_lat[0] + radius, lon_lat[1] - radius],
[lon_lat[0] - radius, lon_lat[1] - radius],
[lon_lat[0] - radius, lon_lat[1] + radius],
[lon_lat[0] + radius, lon_lat[1] + radius]]
)
if area_tyoe == 'roundness':
return Point(lon_lat).buffer(radius)
def roundness(points):
'''
最小外圆形:
首先计算最小外凸多边形，然后计算质心点
然后计算质心点与最小外凸多边形各个顶点
的最远距离，即为对应的最小外围圆范围
:param points: shapely.geometry.multipoint.MultiPoint
:return: center_point:list, roundness:shapely.geometry.polygon.Polygon
'''
polys = points.convex_hull
border = polygon2list(polys) # 最小外凸多边形
center_point = polys.centroid.to_wkt().replace('POINT (','').replace(')','').split(' ') # 质心点
center_point = [round(float(dot),6)for dot in center_point]
buffer_length = [ Point(center_point).distance(Point(point)) for point in border]
buffer_length = round(max(buffer_length),6)
return center_point,Point(center_point).buffer(buffer_length)
def drawPic(points,polys,file_path):
'''
画图，作为开发中参考画出的图形的实际效果
:param points: 点
:param polys: 图形边界
:param file_path: 存储文件位置
:return: 无
'''
from descartes.patch import PolygonPatch
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
fig = plt.figure() # 设定图片大小
ax = fig.add_subplot()
patch = PolygonPatch(polys, facecolor='#6495ED', edgecolor='#6495ED', alpha=0.5, zorder=2)
ax.add_patch(patch)
for point in points:
ax.plot(point.x, point.y, 'o', color='GRAY')
# plt.show()
plt.savefig(file_path)
plt.close()
# **************************边界计算****************************
# 计算每个聚落的最小外凸多边形，最小外矩形，最小外圆形
splace_data_border = pd.DataFrame()
for label in splace_data_lables:
# label = 60
temp = splace_data[splace_data['label'] == label]
points = MultiPoint(temp[['longitude', 'latitude']].values)
# 存在数据质量问题，需要判断全列是否都是一个数
if len(temp['longitude'].unique()) * len(temp['latitude'].unique()) == 1:
# 当数据全一样时，直接以数据点作为原点，100米作为半径
lon_lat = [list(temp['longitude'].unique())[0], list(temp['latitude'].unique())[0]]
polys = createArea(lon_lat, area_tyoe='roundness')
border = polygon2list(polys)
center_point = lon_lat
# 存在数据质量问题，需要判断全列是否都是成一条直线
elif 'LINESTRING' in points.convex_hull.to_wkt():
# 当数据是一条直线时，取其中点作为原点，其中点到一端得距离为半径
temp_str = points.convex_hull.to_wkt().replace('LINESTRING ','').replace('(','').replace(')','')
temp_str = temp_str.split(',')
temp_str = [i.strip().split(' ') for i in temp_str]
lon_lat = [(float(temp_str[0][0]) + float(temp_str[1][0]))/2,
(float(temp_str[0][1]) + float(temp_str[1][1]))/2]
radius = round(Point(lon_lat).distance(Point([float(a) for a in temp_str[0]])),6)
# 20200629添加对于半径过小导致精度不满足，计算结果为0情况修正
if radius == 0:
radius = 0.000009
polys = createArea(lon_lat, area_tyoe='roundness',radius = radius)
border = polygon2list(polys)
center_point = lon_lat
else:
# 画图可以注释掉
polys = points.convex_hull
# drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_convex_hull.jpg')
border_3 = polygon2list(polys) # 最小外凸多边形
polys = points.envelope
# drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_envelope.jpg')
border_1 = polygon2list(polys) # 最小外矩形(边与坐标轴平行)
polys = points.minimum_rotated_rectangle
# drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_minimum_rotated_rectangle.jpg')
border_2 = polygon2list(polys) # 最小外矩形(边不限于与坐标轴平行)
polys = roundness(points)
# drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_roundness.jpg')
border = polygon2list(polys) # 最小外圆形
4145
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
