4.1心得

本文介绍了如何使用PCA方法对地理空间数据中的轮廓线转折点进行降维,然后通过动态规划求解不同建筑物之间的最短路径,并可视化结果。
摘要由CSDN通过智能技术生成

得到沿线生成点,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()

得到

至此,短暂结束

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值