请给以下代码做好排版:
import numpy as np
import rasterio
import matplotlib.pyplot as plt
import pandas as pd
from math import sqrt, atan2, degrees, radians, sin, cos , arccos
import warnings
warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning)
# 1. 读取地形数据
def read_terrain_data(file_path):
with rasterio.open(file_path) as src:
elevation_data = src.read(1)
return elevation_data
# 2. 获取栅格高程
def get_elevation(x, y, elevation_map):
rows, cols = elevation_map.shape
col_idx = x
row_idx = rows - 1 - y
[r,c]=(np.array([row_idx,col_idx])).astype(np.int16)
if 0 <= col_idx < cols and 0 <= row_idx < rows:
return elevation_map[r,c]
else:
return None
# 3. 计算坡度(中心差分法)
def calculate_slope(x, y, elevation_map,flag):
h = get_elevation(x, y, elevation_map)
if h is None:
return None
h_d = get_elevation(x - 1, y, elevation_map)
h_f = get_elevation(x + 1, y, elevation_map)
h_b = get_elevation(x, y + 1, elevation_map)
h_h = get_elevation(x, y - 1, elevation_map)
h_a = get_elevation(x - 1, y+1, elevation_map)
h_c = get_elevation(x + 1, y+1, elevation_map)
h_g = get_elevation(x - 1, y-1, elevation_map)
h_i = get_elevation(x + 1, y-1, elevation_map)
if None in [h_a,h_b,h_c,h_d,h_f,h_g,h_h,h_i]:
return None
dz_dx = (h_c+2*h_f+h_i-h_a-2*h_d-h_g)/8
dz_dy = (h_a+2*h_b+h_c-h_g-2*h_h-h_i)/8
slope_rad = atan2(sqrt(dz_dx**2 + dz_dy**2), 1) #S
slope_deg = degrees(slope_rad)
slope_way = 270 + degrees(dz_dx/dz_dy) - 90*(dz_dy/abs(dz_dy))
way_rad = radians(slope_way) #A
[n1,n2,n3] = [sin(slope_rad)*sin(way_rad),sin(slope_rad)*cos(way_rad),cos(slope_rad)]
if flag== 0:
return slope_deg
else:
return [n1,n2,n3]
# 4. 读取路径数据
def read_path_data(excel_path):
df = pd.read_excel(excel_path, header=None)
da1 = pd.read_excel("D:\model\area1.xlsx",sheet_name=0,header=1)
path1 =[]
da2 = pd.read_excel("D:\model\area1.xlsx",sheet_name=1,header=1)
path2 =[]
for _, row in da1.iterrows():
x, y = int(row[2]), int(row[3])
path1.append((x, y))
for _, row in da2.iterrows():
x, y = int(row[2]), int(row[3])
path2.append((x, y))
path_data = []
for _, row in df.iterrows():
x, y, heading = int(row[1]), int(row[2]), int(row[3])
if 8300<=x<=9800 or 3501<=x<=4500:
if (x,y)in path1 or (x,y)in path2:
path_data.append((x, y, heading,1))
else:
path_data.append((x, y, heading,0))
else:
path_data.append((x, y, heading,0))
return path_data
# 5. 计算路径指标
def calculate_path_metrics(path_data, elevation_map):
total_distance = 0
all_time = 0
dangerous_time =0
esplion =0
angle_changes = []
speeding = [] #速度
length = [] #里程
timing = [0] #时间
heigh = [] #高程
sloping = [] #坡度
for i in range(1, len(path_data)):
x1, y1, h1 = path_data[i - 1]
x2, y2, h2, flag = path_data[i]
heigh.append(h2)
dx = abs(x2 - x1)
dy = abs(y2 - y1)
dl = dx +dy
dheading = abs(h1 -h2)
if dl==1:
if dheading == 0:
w = 1
if dheading == 45:
w = 1.5
if dheading == 90:
w = 2
else:
if dheading == 0:
w = sqrt(2)
if dheading == 45:
w = sqrt(2) + 0.5
if dheading == 90:
w = sqrt(2) + 1
distance = 5*w
length.append(distance) #里程表
total_distance += distance
slope1 = calculate_slope(x1, y1, elevation_map)
slope2 = calculate_slope(x2, y2, elevation_map)
sloping.append(slope2) #坡度表
delta_heading = (h2 - h1 + 180) % 360 - 180
angle_changes.append(abs(delta_heading))
if slope2 is None or slope2 > 30:
speed = 0
else:
speed = 30 if slope2 < 10 else (20 if slope2 < 20 else 10)
speeding.append(speed) #速度表
time = distance / (speed / 3.6) if speed > 0 else 0
timing.append(time) #时间表
all_time += time #时效性
if flag==1:
dangerous_time+=time #安全性
dslope = (slope1 + slope2)/2
n1,n2,n3=calculate_slope(x1, y1, elevation_map,1)
n4,n5,n6=calculate_slope(x2, y2, elevation_map,1)
delta = arccos((n1*n4+n2*n5+n3*n6)/(sqrt(n1**2+n2**2+n3**2)*sqrt(n4**2+n5**2+n6**2)))
esplion += delta * dslope #平稳性
return {
"里程(m)": total_distance,
"时效性(秒)": all_time,
"平稳性": esplion,
"安全性": dangerous_time
}
# 6. 绘制曲线
def plot_curves(a, elevation_map):
# 计算各项指标(保持不变)
distances = [0]
elevations = [get_elevation(path_data[0][0], path_data[0][1], elevation_map)]
slopes = []
speeds = []
angle_changes = [0]
for i in range(1, len(path_data)):
x1, y1, h1 = path_data[i - 1]
x2, y2, h2 = path_data[i]
dx = (x2 - x1) * 5
dy = (y2 - y1) * 5
distance = distances[-1] + sqrt(dx**2 + dy**2)
distances.append(distance)
elevation = get_elevation(x2, y2, elevation_map)
elevations.append(elevation)
slope = calculate_slope(x2, y2, elevation_map)
slopes.append(slope if slope is not None else 0)
delta_heading = (h2 - h1 + 180) % 360 - 180
angle_changes.append(abs(delta_heading))
if slope is None :
speed = 0
else:
base_speed = 30 if slope < 10 else (20 if slope < 20 else 10)
speed = max(5, base_speed - abs(delta_heading) / 10)
speeds.append(speed)
# 转换为numpy数组
distances = np.array(distances)
speeds = np.array(speeds)
angle_changes = np.array(angle_changes)
# 优化3:避免降采样,改用抗锯齿绘制(关键改进)
fig, axes = plt.subplots(4, 1, figsize=(12, 16), dpi=120) # 提高DPI使线条更细
# 1. 高程-里程曲线(保持不变)
axes[0].plot(distances, elevations, 'b-', linewidth=1, alpha=0.8)
axes[0].set_title("Elevation-mileage curve")
axes[0].set_ylabel("Elevation(m)")
# 2. 坡度-里程曲线(保持不变)
axes[1].plot(distances[1:], slopes, 'r-', linewidth=1, alpha=0.8)
axes[1].set_title("Slope-mileage curve")
axes[1].set_ylabel("Slope(°)")
# 3. 速度-距离曲线(优化为细曲线)
axes[2].plot(distances[1:], speeds, 'g-', linewidth=1.2, alpha=0.9,
rasterized=True) # 启用抗锯齿
axes[2].set_title("Speed-distance curve")
axes[2].set_ylabel("Speed(km/h)")
axes[2].set_ylim(0, 40)
axes[2].grid(True, linestyle='--', alpha=0.5)
# 4. 转向角-里程曲线(优化为细曲线)
axes[3].plot(distances, angle_changes, 'm-', linewidth=1.2, alpha=0.9,
rasterized=True) # 启用抗锯齿
axes[3].set_title("Steering angle-mileage curve")
axes[3].set_ylabel("Steering angle(°)")
axes[3].set_xlabel("Distance (m)")
axes[3].set_ylim(0, 40)
axes[3].grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
# 主程序
if __name__ == "__main__":
elevation_map = read_terrain_data("D:/model/map.tif")
path_data = read_path_data("D:/model/pathy.xlsx")
metrics = calculate_path_metrics(path_data, elevation_map)
print("路径指标:")
for key, value in metrics.items():
print(f"{key}: {value:.2f}")
plot_curves(path_data, elevation_map)