点赞发Nature
关注中Science
上一版是直接用了起始位置 + 运输量，代码在这：
基于距离矩阵的OD图
但是我们还希望加入交通运输轨迹，代码如下：
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import pickle
import seaborn as snspath = r"E:\tencent files\chrome Download\Research\LCA\LCA Coal_fired power plant phaseout\data"
map_path = r"E:\tencent files\chrome Download\QGIS 3.12\中国基础数据\中国行政区.shp"trans_2014 = r"coal_transport2014"
geo_data_file = r"Distance_matrix"
transport_matrix = r"coal_trasport2014.xlsx"
mines = "node_mine_pair.xlsx"
power_plants = "node_power_pair.xlsx"
file = r"Results.xlsx"
transport_route_path = os.path.join(path, trans_2014, geo_data_file)
transport_quantity_path = os.path.join(path, trans_2014, transport_matrix)def read_gdf_file(transport_route_path, i):#读取运输轨迹的pickle文件os.chdir(transport_route_path)gdf_file_list = os.listdir(transport_route_path)with open(gdf_file_list[i], "rb") as file:transport_gdf = pickle.load(file)return transport_gdfdef read_distance_matrix(transport_quantity_path):#读取运输量的二维矩阵图trans_mat = pd.read_excel(transport_quantity_path, index_col=0).T.stack()return trans_matdef gdf_concat(transport_route_path, transport_quantity_path, i):transport_quantity_df = read_distance_matrix(transport_quantity_path)gdf_file_list = os.listdir(transport_route_path)transport_quan_pp = transport_quantity_df[int(gdf_file_list[i][:6])]transport_route_gdf = read_gdf_file(transport_route_path, i).set_index(transport_quan_pp.index)transport_df = pd.concat([transport_route_gdf, transport_quan_pp], axis=1)transport_df_route = transport_df[transport_df.iloc[:, -1] != 0]return transport_df_routedef gdf_total(transport_route_path, transport_quantity_path, pp_numbers):#读取每个起止位置的运输量ls = []for i in range(pp_numbers):transport_df_route = gdf_concat(transport_route_path, transport_quantity_path, i)ls.append(transport_df_route)gdf_total = pd.concat(ls)gdf_total.crs = "epsg:32643"gdf_total = gdf_total.to_crs("epsg:4326")return gdf_totaldef CHN_map(map_path):# 读取中国行政区域地图CHN_adm = gpd.read_file(map_path)CHN_adm.to_crs("epsg:4326")return CHN_admdef map_draw(transport_route_path, transport_quantity_path, pp_numbers, map_path=map_path
):fig, ax = plt.subplots(figsize=(16, 16))plt.axis("off")plt.setp(ax.get_yticklabels(), visible=False)plt.setp(ax.get_xticklabels(), visible=False)plt.xlim((90, 140))plt.ylim((30, 50))world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))world.to_crs(epsg=4326)world.plot(ax=ax, color="black", edgecolor="black")CHN_adm = CHN_map(map_path)CHN_adm.plot(ax=ax, color="black", edgecolor="black", linewidth=3)CHN_adm.plot(ax=ax, color="black", edgecolor="white", linewidth=2)transport_gdf = gdf_total(transport_route_path, transport_quantity_path, pp_numbers)transport_gdf.plot(ax=ax)map_draw(transport_route_path, transport_quantity_path, len(os.listdir(transport_route_path))
)
效果图如下：
有不明白的可以私信联系我交流。