得到沿线生成点,buildings1和2的点的个数分别为60和35
接着运行1降维.py
# 输入为轮廓点的坐标,经过PCA方法将二维坐标降维到一维
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
from shapely.geometry import Polygon,box,LineString
import math
from osgeo import ogr
from matplotlib.font_manager import FontProperties
import shapely
import numpy as np
from sklearn.decomposition import PCA
# 读取轮廓线转折点之后的shp文件
point_shp = r'C:\Users\m1959\Desktop\buildings2_all_data\8outline_to_points\output_lines_GeneratePointsAlongLines.shp'
point_gdf = gpd.read_file(point_shp)
# 提取点的坐标数据
points = np.column_stack((point_gdf.geometry.x, point_gdf.geometry.y))
# 使用PCA进行降维
pca = PCA(n_components=1)
pca.fit(points)
# 获取降维后的数据
transformed_points = pca.transform(points)
# 将降维后的数据添加到原始geopandas dataframe
point_gdf['pca_1d'] = transformed_points
# 导出带有PCA数据的shp文件
point_gdf.to_file(r'C:\Users\m1959\Desktop\buildings2_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA.shp')
pd.set_option('display.max_rows', None) # 显示全部行
pd.set_option('display.max_columns', None) # 显示全部列
# 可以通过以下方式查看数据
print(point_gdf)
接着新建ID列,经过一定规则,将拥有pca属性的文件中的排序按照ID升序排列
# 绘制以ID列作为横坐标,pca_1d作为纵坐标的图像并进行可视化
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# 读取shp文件
gdf1 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings1_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA.shp')
# 读取shp文件
gdf2 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings2_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA.shp')
# # 读取shp文件
# gdf3 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings3_all_data\4buildings3_outline\buildings3_Agg_FeatureVertic_withPCA.shp')
# print(gdf1)
# 按ID列升序排列
gdf_1 = gdf1.sort_values(by='ID')
gdf_2 = gdf2.sort_values(by='ID')
# gdf_3 = gdf3.sort_values(by='ID')
pd.set_option('display.max_rows', None) # 显示全部行
pd.set_option('display.max_columns', None) # 显示全部列
# 当ID列第一行中pca_1d大于0时,将整列pca_1d的值取其相反数
if gdf_1.loc[gdf_1['ID'] == 0, 'pca_1d'].values[0] > 0:
gdf_1['pca_1d'] = -1 * gdf_1['pca_1d']
# 当ID列第一行中pca_1d大于0时,将整列pca_1d的值取其相反数
if gdf_2.loc[gdf_2['ID'] == 0, 'pca_1d'].values[0] > 0:
gdf_2['pca_1d'] = -1 * gdf_2['pca_1d']
# # 当ID列第一行中pca_1d大于0时,将整列pca_1d的值取其相反数
# if gdf_3.loc[gdf_3['ID'] == 0, 'pca_1d'].values[0] > 0:
# gdf_3['pca_1d'] = -1 * gdf_3['pca_1d']
# 创建散点图
fig, ax = plt.subplots()
ax.scatter(gdf_1['ID'], gdf_1['pca_1d'], label='buildings1')
ax.scatter(gdf_2['ID'], gdf_2['pca_1d'], label='buildings2')
# ax.scatter(gdf_3['ID'], gdf_3['pca_1d'], label='buildings3')
# 将点按顺序连接起来,不闭合底部
ax.plot(gdf_1['ID'], gdf_1['pca_1d'], label='buildings', marker='o')
# 将点按顺序连接起来,不闭合底部
ax.plot(gdf_2['ID'], gdf_2['pca_1d'], label='buildings', marker='o')
# 将点按顺序连接起来,不闭合底部
# ax.plot(gdf_3['ID'], gdf_3['pca_1d'], label='buildings', marker='o')
# 添加图例
plt.legend()
# 显示图形
plt.show()
# 按ID列升序排列的属性表导出到另一个shp文件中
gdf_1.sort_values(by='ID').to_file(r'C:\Users\m1959\Desktop\buildings1_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA_ID.shp')
gdf_2.sort_values(by='ID').to_file(r'C:\Users\m1959\Desktop\buildings2_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA_ID.shp')
# gdf_3.sort_values(by='ID').to_file(r'C:\Users\m1959\Desktop\buildings3_all_data\4buildings3_outline\buildings3_Agg_FeatureVertic_withPCA_ID.shp')
运行得到以下图像
接着运行3最短路径.py
# 距离矩阵生成累计距离矩阵cumulative_distance_matrix,接下来使用动态规划方法求出最短路径并显示路径
# arr1是纵轴,arr2是横轴
# 经过实验验证合理
import numpy as np
from scipy.spatial import distance
import networkx as nx
import matplotlib.pyplot as plt
import geopandas as gpd
# 读取shp文件
gdf1 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings1_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA_ID.shp')
# 读取shp文件
gdf2 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings2_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA_ID.shp')
# # 读取shp文件
# gdf3 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings3_all_data\4buildings3_outline\buildings3_Agg_FeatureVertic_withPCA_ID.shp')
# 按ID列升序排列
gdf_1 = gdf1.sort_values(by='ID')
gdf_2 = gdf2.sort_values(by='ID')
# gdf_3 = gdf3.sort_values(by='ID')
# 将pca_1d属性的值导出为一维序列
arr1 = np.array(gdf1['pca_1d'])
arr2 = np.array(gdf2['pca_1d'])
# arr3 = np.array(gdf3['pca_1d'])
print("序列1为",arr1)
print("序列2为",arr2)
# print("序列3为",arr3)
# 构造距离矩阵,arr1是纵轴,arr2是横轴
distance_matrix = distance.cdist(arr1.reshape(-1, 1), arr2.reshape(-1, 1), metric='euclidean')
# 可视化距离矩阵为热图
plt.imshow(distance_matrix, cmap='viridis', interpolation='nearest')
plt.colorbar()
plt.show()
# 计算累计距离矩阵
cumulative_distance_matrix = np.zeros_like(distance_matrix)
cumulative_distance_matrix[0, 0] = distance_matrix[0, 0] # 第一个元素与自身的距离为0
for i in range(1, distance_matrix.shape[0]):
cumulative_distance_matrix[i, 0] = cumulative_distance_matrix[i - 1, 0] + distance_matrix[i, 0]
for j in range(1, distance_matrix.shape[1]):
cumulative_distance_matrix[0, j] = cumulative_distance_matrix[0, j - 1] + distance_matrix[0, j]
for i in range(1, distance_matrix.shape[0]):
for j in range(1, distance_matrix.shape[1]):
cumulative_distance_matrix[i, j] = distance_matrix[i, j] + min(cumulative_distance_matrix[i - 1, j], cumulative_distance_matrix[i, j - 1], cumulative_distance_matrix[i - 1, j - 1])
print("累积距离矩阵为",cumulative_distance_matrix)
# 定义斜向移动的权重
diagonal_weight = 1 # 假设斜向移动的权重为1
# 获取累积距离矩阵的维度
m, n = cumulative_distance_matrix.shape
print(m)
print(n)
# 创建路径矩阵,并初始化第一行和第一列
path_matrix = np.zeros((m, n))
path_matrix[0, 0] = cumulative_distance_matrix[0, 0]
for i in range(1, m):
path_matrix[i, 0] = path_matrix[i-1, 0] + cumulative_distance_matrix[i, 0]
for j in range(1, n):
path_matrix[0, j] = path_matrix[0, j-1] + cumulative_distance_matrix[0, j]
# 动态规划求解最短路径
for i in range(1, m):
for j in range(1, n):
# 对于每个位置(i, j),更新路径矩阵的值
# 考虑三种不同的移动方向
move_right = path_matrix[i, j-1] + cumulative_distance_matrix[i, j]
move_down = path_matrix[i-1, j] + cumulative_distance_matrix[i, j]
move_diagonal = path_matrix[i-1, j-1] + cumulative_distance_matrix[i, j]*diagonal_weight
path_matrix[i, j] = min(move_right, move_down, move_diagonal)
# 输出起点到终点的最短路径
print("最短路径距离为:", path_matrix[m-1, n-1])
# 反向推导出最短路径
i, j = m-1, n-1
shortest_path = [(i, j)]
while i > 0 or j > 0:
move_right = path_matrix[i, j-1]
move_down = path_matrix[i-1, j]
move_diagonal = path_matrix[i-1, j-1]
min_move = min(move_right, move_down, move_diagonal)
if min_move == move_right:
j -= 1
elif min_move == move_down:
i -= 1
else:
i -= 1
j -= 1
shortest_path.append((i, j))
shortest_path.reverse()
print("最短路径为:", shortest_path)
得到
以及
最短路径距离为: 129502.67837982874
最短路径为: [(0, 0), (0, 1), (0, 2), (0, 3), (1, 4), (2, 5), (2, 6), (3, 7), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12), (9, 13), (10, 14), (11, 15), (12, 16), (13, 17), (14, 18), (15, 19), (16, 20), (17, 21), (18, 22), (19, 23), (20, 24), (21, 25), (22, 26), (23, 27), (24, 27), (25, 27), (26, 28), (27, 28), (28, 28), (29, 28), (30, 28), (31, 28), (32, 28), (33, 28), (34, 28), (35, 28), (36, 28), (37, 28), (38, 28), (39, 28), (40, 28), (41, 28), (42, 28), (43, 28), (44, 28), (45, 28), (46, 28), (47, 28), (48, 28), (49, 28), (50, 28), (51, 28), (52, 28), (53, 28), (54, 29), (55, 30), (56, 31), (57, 32), (58, 33), (59, 34)]
接着运行4填充矩阵.py
# 得到元素个数相等的两个经过匹配校正后的序列
import numpy as np
from scipy.spatial import distance
import networkx as nx
import matplotlib.pyplot as plt
import geopandas as gpd
# 读取shp文件
gdf1 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings1_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA_ID.shp')
# 读取shp文件
gdf2 = gpd.read_file(r'C:\Users\m1959\Desktop\buildings2_all_data\8outline_to_points\outline_GeneratePointsAlongLines_withPCA_ID.shp')
# 按ID列升序排列
gdf_1 = gdf1.sort_values(by='ID')
gdf_2 = gdf2.sort_values(by='ID')
# print(gdf_1)
# 将pca_1d属性的值导出为一维序列
arr1 = np.array(gdf1['pca_1d'])
arr2 = np.array(gdf2['pca_1d'])
print("序列1为",arr1)# 60个元素
print("序列2为",arr2)# 35个元素
# 最短路径一共64步
# 创建一个全零的 60*64 矩阵
matrix1 = np.zeros((60, 64))
coordinates = [(0, 0), (0, 1), (0, 2), (1, 3), (2, 4), (2, 5), (2, 6), (3, 7), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12), (9, 13), (10, 14), (11, 15), (12, 16), (13, 17), (14, 18), (15, 19), (16, 20), (17, 21), (18, 22), (19, 23), (20, 24), (21, 25), (22, 26), (23, 26), (24, 26), (25, 26), (26, 27), (27, 27), (28, 27), (29, 27), (30, 27), (31, 27), (32, 27), (33, 27), (34, 27), (35, 27), (36, 27), (37, 27), (38, 27), (39, 27), (40, 27), (41, 27), (42, 27), (43, 27), (44, 27), (45, 27), (46, 27), (47, 27), (48, 27), (49, 27), (50, 27), (51, 27), (52, 27), (53, 28), (54, 29), (55, 30), (56, 31), (57, 32), (58, 33), (59, 34)]
# 用于存储结果的字典
result1 = {}
# 遍历并重新编号
for i, (x, y) in enumerate(coordinates):
result1[x, i] = 1 # 仅保留重新编号后的结果
# 输出结果
print(result1)
# 使用结果字典中的键将矩阵中对应位置的值设置为1
for key in result1:
matrix1[key] = 1
# 输出结果
print("规整矩阵1为",matrix1)
# 创建一个 64*35 全零矩阵
matrix2 = np.zeros((64, 35))
# 将coordinates中的行重新编号为0-63
for i, (x, y) in enumerate(coordinates):
matrix2[i, y] = 1
# 设置numpy打印选项,以显示所有元素
np.set_printoptions(threshold=np.inf)
#求转置矩阵
transposed_matrix = matrix2.T
print(transposed_matrix)
# 输出结果
print("规整矩阵2为",transposed_matrix)
print("矩阵的行数:", matrix2.shape[0])
# 计算 arr1 和 matrix1 的点积
dot_product1 = np.dot(arr1, matrix1)
dot_product2 = np.dot(arr2, transposed_matrix)
# 打印点积
print("1的点积为:",dot_product1)
print("2的点积为:",dot_product2)
# 创建图表
plt.plot(dot_product1, label="dot_product1")
plt.plot(dot_product2, label="dot_product2")
# 添加标签
plt.xlabel("Index")
plt.ylabel("Value")
plt.title("Dot Product")
# 添加图例
plt.legend()
# 显示图表
plt.show()
得到
至此,短暂结束