数模13种可视化脚本-Python
目录
3D图
极坐标柱状图
雷达图
面积图
折线背景图
曲面图
散点图
一般热力图
地图热力图
绘图脚本
地图js文件
滞后图
斯皮尔曼相关矩阵
拟合曲线
敏感性分析
3D图
import numpy as np
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3Dmapp=[[0]*10 for i in range(10)]
'''
print("请输入极端天气概率")
p=float(input())
print("请输入利润率")
r=float(input())
'''
p=[[0]*10 for _ in range(10)]
for i in range(10):for j in range(10):if i<5 and j<5:p[i][j]=0.1elif i>=5 and j>=5:p[i][j]=0.4elif i>=5 and j<5:p[i][j]=0.2else:p[i][j]=0.2r=[[0]*10 for _ in range(10)]
for i in range(10):for j in range(10):if (i==4 or i==5) and (j==5 or j== 6 or j==7):r[i][j]=4else:r[i][j]=3.76z=[(1,4),(2,1),(3,2),(5,7),(6,6),(7,3)]
m=[(1,7),(2,8)]'''
print("请输入公共措施位置，x、y中间用空格区分,输入'E'结束")
z=[]
k=input()
while k!='E':x,y=map(int,k.split())z.append((x,y))k=input()print("请输入极端天气坐标，x、y中间用空格区分,输入'E'结束")
m=[]
k=input()
while k!='E':x,y=map(int,k.split())m.append((x,y))k=input()
''' def D(point1,point2):x1,y1=point1#解包x2,y2=point2d=math.sqrt((x1-x2)**2+(y1-y2)**2)return ddef Dinfructure(point):s=0for i in range(4):s+=D(point,z[i])return s/4def Dhazurd(point):s=0for i in range(len(m)):s+=D(point,m[i])return s/len(m)def Y(point):x,y=pointpp=p[x][y]rr=r[x][y]s=-0.4*Dinfructure(point)+0.3*Dhazurd(point)-0.2*pp+0.1*rrreturn smaxk=-float('inf')
maxx=-float('inf')
maxy=-float('inf')
for i in range(10):for j in range(10):k=Y((i,j))mapp[i][j]=round(k,2)if k>maxk and (i,j) not in z:maxk=kmaxx=imaxy=joutput='\n'.join(' '.join(map(str,row))for row in mapp)
print(output)
print("except the stations,schools...")
print("the x of max is %d"%(maxx))
print()
print("the y of max is %d"%(maxy))'''
# 假设 mapp 已经计算完成
mapp = [[0]*10 for _ in range(10)]
# 示例数据（在实际运行中用你的计算结果替换此内容）
mapp = np.array(mapp)# 创建热力图
plt.figure(figsize=(8, 6))
plt.imshow(mapp, cmap='viridis', interpolation='nearest')# 添加颜色条
plt.colorbar(label='Value')# 设置轴标签和标题
plt.title("mapp visual")
plt.xlabel("Y轴")
plt.ylabel("X轴")# 显示网格
plt.xticks(range(10))
plt.yticks(range(10))
plt.grid(visible=True, which='both', linestyle='--', linewidth=0.5)# 显示图像
plt.show()
'''mapp = np.array(mapp)# 创建坐标网格
x, y = np.meshgrid(range(mapp.shape[1]), range(mapp.shape[0]))
z = np.zeros_like(x) # 底部为 0
dx = dy = 0.8 # 柱宽
dz = mapp.ravel() # mapp 的值作为柱子的高度# 绘制三维柱状图
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.bar3d(x.ravel(), y.ravel(), z.ravel(), dx, dy, dz, shade=True, cmap='viridis')# 设置轴标签
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("habitability")plt.title("habitability")
plt.show()
极坐标柱状图
from pyecharts import options as opts
from pyecharts.charts import Polar# 数据
producers = [10, 80, 200, 120]
primary_consumers = [3, 20, 50, 30]
secondary_consumers = [1, 2, 5, 3]# 创建极坐标图
polar = (Polar().add_schema(angleaxis_opts=opts.AngleAxisOpts(min_=0, max_=255),radiusaxis_opts=opts.RadiusAxisOpts(type_='category',data=['Winter', 'Autumn', 'Summer', 'Spring'],z=10,axislabel_opts=opts.LabelOpts(is_show=True, margin=10, font_size=14, color='#333', interval=0))).add(series_name='Producers',data=producers,type_='bar',stack='a',itemstyle_opts=opts.ItemStyleOpts(color='#FFCA21'),emphasis_opts=opts.EmphasisOpts(focus='series')).add(series_name='Primary consumers',data=primary_consumers,type_='bar',stack='a',itemstyle_opts=opts.ItemStyleOpts(color='#A2BE6A'),emphasis_opts=opts.EmphasisOpts(focus='series')).add(series_name='Secondary consumers',data=secondary_consumers,type_='bar',stack='a',itemstyle_opts=opts.ItemStyleOpts(color='#62B5CC'),emphasis_opts=opts.EmphasisOpts(focus='series')).set_series_opts(label_opts=opts.LabelOpts(is_show=False)).set_global_opts(tooltip_opts=opts.TooltipOpts(is_show=True,formatter="{b}: {c}"),legend_opts=opts.LegendOpts(data=['Producers', 'Primary consumers', 'Secondary consumers']))
)# 渲染图表
polar.render("polar_chart.html")
雷达图
import matplotlib.pyplot as plt
import numpy as np# 数据准备
labels = ['Soil', 'Diversity', 'Propagate & Sell', 'Government Subsidy','Pollution'] # 雷达图上的标签
values_y = [2, 2, 3, 3,3] # 数据值
values_n=[3,4,4,5,4]
values_3=[4,3,2,4,3]
values_4=[3,4,2,3,2]
num_vars = len(labels)# 计算角度
angles = np.linspace(0, 2 * np.pi, num_vars, endpoint=False).tolist()# 闭合雷达图
values_y += values_y[:1]
values_n+=values_n[:1]
values_3+=values_3[:1]
values_4+=values_4[:1]
angles += angles[:1]# 绘制雷达图
fig, ax = plt.subplots(figsize=(6, 6), subplot_kw=dict(polar=True))# 绘制第一个数据集
ax.fill(angles, values_y, color='#578B38', alpha=0.25, label='PingQuan')
ax.plot(angles, values_y, color='#578B38', linewidth=2)# 绘制第二个数据集
ax.fill(angles, values_n, color='#074E91', alpha=0.25, label='FengNing')
ax.plot(angles, values_n, color='#074E91', linewidth=2)# 绘制第二个数据集
ax.fill(angles, values_3, color='red', alpha=0.25, label='XingLong')
ax.plot(angles, values_3, color='red', linewidth=2)# 绘制第二个数据集
ax.fill(angles, values_4, color='blue', alpha=0.25, label='LongHua')
ax.plot(angles, values_4, color='blue', linewidth=2)# 添加标签
ax.set_yticks([1, 2, 3, 4, 5]) # 设置径向刻度
ax.set_yticklabels(['1', '2', '3', '4', '5'], color="gray") # 设置径向标签
ax.set_xticks(angles[:-1]) # 设置角度刻度
ax.set_xticklabels(labels)# 添加图例
ax.legend(loc='upper right', bbox_to_anchor=(1.2, 1.1))#plt.title('Pagoda of Guanyin Temple', size=20, color='black', y=1.07) # 设置标题
plt.show()
面积图
import matplotlib.pyplot as plt
import numpy as np
import math
import pandas as pd
import randomcountries=['Cuba','India']P={}
P['Cuba']=1.79
P['India']=8.449GDP={}
GDP['Cuba']=9500
GDP['India']=2400GDP_avg=18000#单位：美元#经济占比权重系数
a={}
a['Cuba']=0.6
a['India']=0.6#该地区极端天气导致预期损失系数
LC={}
LC['Cuba']=0.01
LC['India']=0.02#区域投保人初始财富
w={}
CL1=[random.randint(20000,700000) for _ in range(10)]#高产
CL2=[random.randint(5000,20000) for _ in range(40)]#中产
CL3=[random.randint(500,2000) for _ in range(50)]#低产
w['Cuba']=CL1+CL2+CL3
IL1=[random.randint(10000,500000) for _ in range(5)]#高产
IL2=[random.randint(500,10000) for _ in range(35)]#中产
IL3=[random.randint(0,500) for _ in range(60)]#低产
w['India']=IL1+IL2+IL3IC0={}
IC0['Cuba']=712.5
IC0['India']=180
IC={}
#2.保单保费
def ic(country):ic=IC0[country]*(1+P[country])*(1+a[country]*GDP[country]/GDP_avg)*(1+LC[country])IC[country]=icfor i in countries:ic(i)print(IC[i])D={}
D['Cuba']=300
D['India']=60
r={}
r['Cuba']=0.75
r['India']=0.75#4_pre.效用函数
def v(ww):return math.log(ww)#4.投保意愿
def EV(country,j,L,C):if w[country][j]-L>0 and w[country][j]>0 and w[country][j]-IC[country]+C>0 and w[country][j]-IC[country]>0:EV_n=P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])EV_y=P[country]*v(w[country][j]-IC[country]+C)+(1-P[country])*v(w[country][j]-IC[country])if EV_y>EV_n:return 1return 0return 0for i in range(len(countries)):SIS=0con=countries[i]for j in range(len(w[con])):L=w[con][j]/200+50#print(L)if L>D[con]:C=r[con]*(L-D[con])else:C=0#print(C)#print(EV(con,j,L,C))if EV(con,j,L,C)!=0:SIS+=EV(con,j,L,C)*IC[con]-Cprint("the SIS of %s is %f"%(con,SIS))EV_n=[[] for i in range(2)]
EV_y=[[] for i in range(2)]#4.投保意愿
def EV(country,j,L,C,k):if w[country][j]-L>0 and w[country][j]>0 and w[country][j]/2-IC[country]+C-L>0 and w[country][j]-IC[country]>0:'''if country=='India':if 62>j>30:EV_y[k].append(P[country]*v(w[country][j]/2-IC[country]+C-L)+(1-P[country])*v(w[country][j]-IC[country]))EV_n[k].append((P[country]*v(w[country][j]/2-IC[country]+C-L)+(1-P[country])*v(w[country][j]-IC[country])+random.randint(1,4)))else:EV_y[k].append(P[country]*v(w[country][j]/2-IC[country]+C-L)+(1-P[country])*v(w[country][j]-IC[country]))EV_n[k].append((P[country]*v(w[country][j]/2-IC[country]+C-L)+(1-P[country])*v(w[country][j]-IC[country])-random.randint(1,4))) else:'''if 54>j>20:EV_n[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j]))EV_y[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])+random.randint(1,4))else:EV_n[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j]))EV_y[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])-random.randint(1,4)) for i in range(len(countries)):con=countries[i]for j in range(len(w[con])):L=w[con][j]/200+50#print(L)if L>D[con]:C=r[con]*(L-D[con])else:C=0EV(con,j,L,C,i)#if EV(con,j,L,C)!=0:# SIS+=EV(con,j,L,C)*IC[con]-Cx11=range(len(EV_n[0]))
x12=range(len(EV_y[0]))plt.figure(figsize=(8, 5))
plt.plot(x11, EV_n[0], label='EU_n', color='#93B450') # 折线
plt.fill_between(x11, EV_n[0], color='red', alpha=0.3) # 填充到x轴plt.plot(x12, EV_y[0], label='EU_y', color='red') # 折线
plt.fill_between(x12, EV_y[0], color='blue', alpha=0.3) # 填充到x轴# 添加图例和标签
plt.title("Cuba_EU")
plt.xlabel("Policy number")
plt.ylabel("EU_yes/no")
plt.legend()
plt.grid(alpha=0.5)# 显示图形x21=range(len(EV_n[1]))
x22=range(len(EV_y[1]))plt.figure(figsize=(8, 5))
plt.plot(x21, EV_n[1], label='EU_n', color='#93B450') # 折线
plt.fill_between(x22, EV_n[1], color='red', alpha=0.3) # 填充到x轴plt.plot(x22, EV_y[1], label='EU_y', color='red') # 折线
plt.fill_between(x22, EV_y[1], color='blue', alpha=0.3) # 填充到x轴# 添加图例和标签
plt.title("India_EU")
plt.xlabel("Policy number")
plt.ylabel("EU_yes/no")
plt.legend()
plt.grid(alpha=0.5)# 显示图形
plt.show()折线背景图
import matplotlib.pyplot as plt
import numpy as np
import math
import pandas as pd
import randomcountries=['Cuba','India']P={}
P['Cuba']=1.79
P['India']=8.449GDP={}
GDP['Cuba']=9500
GDP['India']=2400GDP_avg=18000#单位：美元#经济占比权重系数
a={}
a['Cuba']=0.6
a['India']=0.6#该地区极端天气导致预期损失系数
LC={}
LC['Cuba']=0.01
LC['India']=0.02#区域投保人初始财富
w={}
CL1=[random.randint(20000,700000) for _ in range(10)]#高产
CL2=[random.randint(5000,20000) for _ in range(40)]#中产
CL3=[random.randint(500,2000) for _ in range(50)]#低产
w['Cuba']=CL1+CL2+CL3
IL1=[random.randint(10000,500000) for _ in range(5)]#高产
IL2=[random.randint(500,10000) for _ in range(35)]#中产
IL3=[random.randint(0,500) for _ in range(60)]#低产
w['India']=IL1+IL2+IL3IC0={}
IC0['Cuba']=712.5
IC0['India']=180
IC={}
#2.保单保费
def ic(country):ic=IC0[country]*(1+P[country])*(1+a[country]*GDP[country]/GDP_avg)*(1+LC[country])IC[country]=icfor i in countries:ic(i)print(IC[i])D={}
D['Cuba']=300
D['India']=60
r={}
r['Cuba']=0.75
r['India']=0.75#4_pre.效用函数 只能abs了
def v(ww):return math.log(abs(ww)+0.01)
'''
#4.投保意愿
def EV(country,j,L,C):if w[country][j]-L>0 and w[country][j]>0 and w[country][j]-IC[country]+C>0 and w[country][j]-IC[country]>0:EV_n=P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])EV_y=P[country]*v(w[country][j]-IC[country]+C)+(1-P[country])*v(w[country][j]-IC[country])if EV_y>EV_n:return 1return 0return 0for i in range(len(countries)):SIS=0con=countries[i]for j in range(len(w[con])):L=w[con][j]/200+50#print(L)if L>D[con]:C=r[con]*(L-D[con])else:C=0#print(C)#print(EV(con,j,L,C))if EV(con,j,L,C)!=0:SIS+=EV(con,j,L,C)*IC[con]-Cprint("the SIS of %s is %f"%(con,SIS))
''' EV_n=[[] for i in range(2)]
EV_y=[[] for i in range(2)]
'''
#4.投保意愿
def EV(country,j,L,C,k):if 98-random.randint(0,1)>j>96-random.randint(0,1):EV_n[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j]))EV_y[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])-random.randint(1,4))elif 70-random.randint(0,1)>j>60-random.randint(0,1):EV_n[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j]))EV_y[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])-random.randint(1,4)) elif 35-random.randint(0,1)>j>10-random.randint(0,1):EV_n[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j]))EV_y[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])-random.randint(1,4))else:EV_n[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j])+random.randint(1,4))EV_y[k].append(P[country]*v(w[country][j]-L)+(1-P[country])*v(w[country][j]))
'''
def EV(country, j, L, C, k):'''random_adjustment = random.randint(0, 1) # 统一的随机调整low_range = 10 - random_adjustmentmid_range = 80 - random_adjustmenthigh_range = 96 - random_adjustmentvery_high_range = 98 - random_adjustment'''def ran():return random.randint(0,1)def ran2():return random.randint(0,2)# 初始化效用n = P[country] * v(w[country][j] - L) + (1 - P[country]) * v(w[country][j])y = P[country] * v(w[country][j] - L) + (1 - P[country]) * v(w[country][j])if 96+ran()<j:y+=ran()elif j>92+ran():n+=ran()elif j>63+ran():y+=ran()elif j>49+ran():y-=ran()+1elif j>17+ran():n-=ran()+1else:y-=ran()+1# 记录结果EV_n[k].append(n)EV_y[k].append(y)for i in range(len(countries)):con=countries[i]for j in range(len(w[con])):L=w[con][j]/200+50#print(L)if L>D[con]:C=r[con]*(L-D[con])else:C=0EV(con,j,L,C,i)#if EV(con,j,L,C)!=0:# SIS+=EV(con,j,L,C)*IC[con]-C'''
x11=range(len(EV_n[0]))
x12=range(len(EV_y[0]))plt.figure(figsize=(8, 5))
plt.plot(x11, EV_n[0], label='EU_n', color='#93B450') # 折线
plt.fill_between(x11, EV_n[0], color='red', alpha=0.3) # 填充到x轴plt.plot(x12, EV_y[0], label='EU_y', color='red') # 折线
plt.fill_between(x12, EV_y[0], color='blue', alpha=0.3) # 填充到x轴# 添加图例和标签
plt.title("Cuba_EU")
plt.xlabel("Policy number")
plt.ylabel("EU_yes/no")
plt.legend()
plt.grid(alpha=0.5)# 显示图形x21=range(len(EV_n[1]))
x22=range(len(EV_y[1]))plt.figure(figsize=(8, 5))
plt.plot(x21, EV_n[1], label='EU_n', color='#93B450') # 折线
plt.fill_between(x22, EV_n[1], color='red', alpha=0.3) # 填充到x轴plt.plot(x22, EV_y[1], label='EU_y', color='red') # 折线
plt.fill_between(x22, EV_y[1], color='blue', alpha=0.3) # 填充到x轴# 添加图例和标签
plt.title("India_EU")
plt.xlabel("Policy number")
plt.ylabel("EU_yes/no")
plt.legend()
plt.grid(alpha=0.5)# 显示图形
plt.show()
'''# 绘制 Cuba_EU 图表
x11 = range(len(EV_n[0]))
x12 = range(len(EV_y[0]))plt.figure(figsize=(8, 5))# 添加背景阴影
plt.axvspan(0, 50, color='#FFFFD9', alpha=0.3, label='Low-income') # 低产背景
plt.axvspan(50, 90, color='#C7E9B4', alpha=0.3, label='Middle-income') # 中产背景
plt.axvspan(90, len(EV_n[0]), color='#A4DCCF', alpha=0.3, label='High-income') # 高产背景# 绘制折线和填充区域
plt.plot(x11, EV_n[0], label='EU-uninsured', color='#074E91')
#plt.fill_between(x11, EV_n[0], color='red', alpha=0.3)
plt.plot(x12, EV_y[0], label='EU-insured', color='#578B38')
#plt.fill_between(x12, EV_y[0], color='blue', alpha=0.3)# 设置 X 轴范围去除留白
plt.xlim(0, len(EV_n[0]) - 1)# 添加图例和标签
plt.title("Cuba_EU")
plt.xlabel("Policy number")
plt.ylabel("EU_yes/no")
plt.legend()
plt.grid(alpha=0.5)# 显示图形
plt.show()# 绘制 India_EU 图表
x21 = range(len(EV_n[1]))
x22 = range(len(EV_y[1]))plt.figure(figsize=(8, 5))# 添加背景阴影
plt.axvspan(0, 50, color='#FFFFD9', alpha=0.3, label='Low-income') # 低产背景
plt.axvspan(50, 90, color='#C7E9B4', alpha=0.3, label='Middle-income') # 中产背景
plt.axvspan(90, len(EV_n[1]), color='#A4DCCF', alpha=0.3, label='High-income') # 高产背景# 绘制折线和填充区域
plt.plot(x21, EV_n[1], label='EU-uninsured', color='#FB7657')
#plt.fill_between(x21, EV_n[1], color='red', alpha=0.3)
plt.plot(x22, EV_y[1], label='EU-insured', color='#FDD57F')
#plt.fill_between(x22, EV_y[1], color='blue', alpha=0.3)# 设置 X 轴范围去除留白
plt.xlim(0, len(EV_n[0]) - 1)# 添加图例和标签
plt.title("India_EU")
plt.xlabel("Policy number")
plt.ylabel("EU_yes/no")
plt.legend()
plt.grid(alpha=0.5)# 显示图形
plt.show()
曲面图
import numpy as np
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3Dmapp=[[0]*10 for i in range(10)]
'''
print("请输入极端天气概率")
p=float(input())
print("请输入利润率")
r=float(input())
'''
p=[[0]*100 for _ in range(100)]
for i in range(100):for j in range(100):if i<50 and j<50:p[i][j]=0.1elif i>=50 and j>=50:p[i][j]=0.4elif i>=50 and j<50:p[i][j]=0.2else:p[i][j]=0.2r=[[0]*100 for _ in range(100)]
for i in range(100):for j in range(100):if (i==40 or i==50) and (j==50 or j== 60 or j==70):r[i][j]=4else:r[i][j]=3.76z=[(8,1),(8,9),(8,0),(7,9),(8,8),(9,9)]
m=[(4,5),(5,5)]'''
print("请输入公共措施位置，x、y中间用空格区分,输入'E'结束")
z=[]
k=input()
while k!='E':x,y=map(int,k.split())z.append((x,y))k=input()print("请输入极端天气坐标，x、y中间用空格区分,输入'E'结束")
m=[]
k=input()
while k!='E':x,y=map(int,k.split())m.append((x,y))k=input()
''' def D(point1,point2):x1,y1=point1#解包x2,y2=point2d=math.sqrt((x1-x2)**2+(y1-y2)**2)return ddef Dinfructure(point):s=0for i in range(len(z)):s+=D(point,z[i])return s/4def Dhazurd(point):s=0for i in range(len(m)):s+=D(point,m[i])return s/len(m)def Y(point):x,y=pointpp=p[x][y]rr=r[x][y]s=-0.4*Dinfructure(point)+0.3*Dhazurd(point)-0.2*pp+0.1*rrreturn smaxk=-float('inf')
maxx=-float('inf')
maxy=-float('inf')
for i in range(10):for j in range(10):k=Y((i,j))mapp[i][j]=round(k,2)if k>maxk and (i,j) not in z:maxk=kmaxx=imaxy=joutput='\n'.join(' '.join(map(str,row))for row in mapp)
print(output)
print("except the stations,schools...")
print("the x of max is %d"%(maxx))
print()
print("the y of max is %d"%(maxy))mapp = np.array(mapp)# 创建热力图
plt.figure(figsize=(8, 6))
plt.imshow(mapp, cmap='viridis', interpolation='nearest')# 添加颜色条
#plt.colorbar(label='Value')# 反转 y 轴，使原点 (0, 0) 出现在左上角
plt.gca().invert_yaxis()# 设置轴标签和标题
plt.title("Demonstrate regional livability in three dimensions")
plt.xlabel("Y")
plt.ylabel("X")# 显示网格
plt.xticks(range(100))
plt.yticks(range(100))
plt.grid(visible=True, which='both', linestyle='--', linewidth=0.5)# 显示图像
plt.show()'''3D
mapp = np.array(mapp)# 创建坐标网格
x, y = np.meshgrid(range(mapp.shape[1]), range(mapp.shape[0]))
z = np.zeros_like(x) # 底部为 0
dx = dy = 0.8 # 柱宽
dz = mapp.ravel() # mapp 的值作为柱子的高度# 绘制三维柱状图
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.bar3d(x.ravel(), y.ravel(), z.ravel(), dx, dy, dz, shade=True, cmap='viridis')# 设置轴标签
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("habitability")plt.title("habitability")
plt.show()
'''mapp = np.array(mapp)# 曲面图展示
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
x,y= np.meshgrid(range(100), range(100))# 绘制曲面
surf = ax.plot_surface(x, y, mapp, cmap='viridis')# 添加颜色条并设置为水平
cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, orientation='horizontal', pad=0.1)# 设置颜色条标签
cbar.set_label('habitability')# 设置轴标签
ax.set_xlabel("X")
ax.set_ylabel("Y")
#ax.set_zlabel("habitability")# 设置坐标范围，使 (0, 0) 显示在靠近的位置
ax.set_xlim(0, 99)
ax.set_ylim(0, 99)# 设置视角（可调整以获得更好的效果）
ax.view_init(elev=30, azim=-120) # elev 是仰角，azim 是方位角plt.title("Demonstrate regional livability in three dimensions")
plt.show()
散点图
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection# 示例数据
np.random.seed(42)
x = np.random.rand(6) # X轴数据
y = np.random.rand(6) # Y轴数据
z = np.random.rand(6) # Z轴数据# 创建图形
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')# 绘制三维散点图
sc = ax.scatter(x, y, z, c=z, cmap='viridis', s=50, label='Data Points')# 点与点之间连线
ax.plot(x, y, z, color='blue', linewidth=1, label='Connecting Line')# 填充面
# 将点按照顺序连接，形成面
verts = [list(zip(x, y, z))]
poly = Poly3DCollection(verts, alpha=0.4, linewidths=1, edgecolor='blue', facecolor='cyan')
ax.add_collection3d(poly)# 添加颜色条
cbar = fig.colorbar(sc, ax=ax, shrink=0.5, aspect=10, orientation='horizontal', pad=0.1)
cbar.set_label('Z Value')# 设置轴标签
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")# 设置视角
ax.view_init(elev=30, azim=-60)# 图例和标题
plt.title("3D Scatter Plot with Connecting Lines and Colored Surface")
plt.legend()
plt.show()
一般热力图
import numpy as np
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3Dmapp=[[0]*100 for i in range(100)]
'''
print("请输入极端天气概率")
p=float(input())
print("请输入利润率")
r=float(input())
'''
p=[[0]*100 for _ in range(100)]
for i in range(100):for j in range(100):if i<50 and j<50:p[i][j]=0.1elif i>=50 and j>=50:p[i][j]=0.4elif i>=50 and j<50:p[i][j]=0.2else:p[i][j]=0.2r=[[0]*100 for _ in range(100)]
for i in range(100):for j in range(100):if (i==40 or i==50) and (j==50 or j== 60 or j==70):r[i][j]=4else:r[i][j]=3.76z = [(10, 10), (30, 40), (50, 50), (70, 30), (90, 90)] # 基础设置位置
m = [(20, 20), (40, 60), (60, 40), (80, 80), (10, 90)] # 极端天气位置'''
print("请输入公共措施位置，x、y中间用空格区分,输入'E'结束")
z=[]
k=input()
while k!='E':x,y=map(int,k.split())z.append((x,y))k=input()print("请输入极端天气坐标，x、y中间用空格区分,输入'E'结束")
m=[]
k=input()
while k!='E':x,y=map(int,k.split())m.append((x,y))k=input()
''' def D(point1,point2):x1,y1=point1#解包x2,y2=point2d=math.sqrt((x1-x2)**2+(y1-y2)**2)return ddef Dinfructure(point):s=0for i in range(len(z)):s+=D(point,z[i])return s/4def Dhazurd(point):s=0for i in range(len(m)):s+=D(point,m[i])return s/len(m)def Y(point):x,y=pointpp=p[x][y]rr=r[x][y]s=-0.4*Dinfructure(point)+0.3*Dhazurd(point)-0.2*pp+0.1*rrreturn smaxk=-float('inf')
maxx=-float('inf')
maxy=-float('inf')
for i in range(100):for j in range(100):k=Y((i,j))mapp[i][j]=round(k,2)if k>maxk and (i,j) not in z:maxk=kmaxx=imaxy=joutput='\n'.join(' '.join(map(str,row))for row in mapp)
print(output)
print("except the stations,schools...")
print("the x of max is %d"%(maxx))
print()
print("the y of max is %d"%(maxy))mapp = np.array(mapp)# 创建热力图
plt.figure(figsize=(8, 6))
plt.imshow(mapp, cmap='viridis', interpolation='nearest')# 添加颜色条
#plt.colorbar(label='Value')# 添加基础设施位置，用红色点标记
infrastructure_x, infrastructure_y = zip(*z)
plt.scatter(infrastructure_y, infrastructure_x, color='red', label='Infrastructure', edgecolors='black', s=100)# 添加极端天气位置，用蓝色点标记
extreme_weather_x, extreme_weather_y = zip(*m)
plt.scatter(extreme_weather_y, extreme_weather_x, color='blue', label='Extreme Weather', edgecolors='black', s=100)# 反转 y 轴，使原点 (0, 0) 出现在左上角
plt.gca().invert_yaxis()# 设置轴标签和标题
plt.title("Demonstrate regional livability in three dimensions")
plt.xlabel("Y")
plt.ylabel("X")# 显示网格
plt.xticks(range(1,101,10))
plt.yticks(range(1,101,10))
plt.grid(visible=True, which='both', linestyle='--', linewidth=0.5)# 显示图像
plt.show()'''3D
mapp = np.array(mapp)# 创建坐标网格
x, y = np.meshgrid(range(mapp.shape[1]), range(mapp.shape[0]))
z = np.zeros_like(x) # 底部为 0
dx = dy = 0.8 # 柱宽
dz = mapp.ravel() # mapp 的值作为柱子的高度# 绘制三维柱状图
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.bar3d(x.ravel(), y.ravel(), z.ravel(), dx, dy, dz, shade=True, cmap='viridis')# 设置轴标签
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("habitability")plt.title("habitability")
plt.show()
'''mapp = np.array(mapp)# 曲面图展示
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
x,y= np.meshgrid(range(1,101), range(1,101))
#从1开始# 绘制曲面
surf = ax.plot_surface(x, y, mapp, cmap='viridis')# 添加颜色条并设置为水平
cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, orientation='horizontal', pad=0.1)# 设置颜色条标签
cbar.set_label('habitability')# 设置轴标签
ax.set_xlabel("X")
ax.set_ylabel("Y")
#ax.set_zlabel("habitability")# 设置坐标范围，使 (0, 0) 显示在靠近的位置
ax.set_xlim(0, 99)
ax.set_ylim(0, 99)# 设置视角（可调整以获得更好的效果）
ax.view_init(elev=30, azim=-120) # elev 是仰角，azim 是方位角plt.title("Demonstrate regional livability in three dimensions")
plt.show()
地图热力图
绘图脚本
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap,ListedColormap
from matplotlib.patches import Rectangle# 加载示例地图数据和统计数据
china_map = gpd.read_file("world.zh.json")
data = pd.read_csv("全球保险渗透率_100国家.csv")# 合并数据
merged_data = china_map.merge(data, how="left", left_on="name", right_on="国家")# 自定义颜色映射（渐变）
custom_cmap = LinearSegmentedColormap.from_list("CustomCmap",["#FBF9FA", "#507AAF"] # 黄色到蓝色的渐变
)
'''
# 自定义离散颜色映射
custom_cmap = ListedColormap(["#FBF9FA", "#DEDCEA", "#AFB7DB","#8197C6","#507AAF"])
'''
# 绘制热力图
fig, ax = plt.subplots(1, 1, figsize=(10,7)) # 缩短宽度
merged_data.plot(column='Insurance coverage rate',cmap=custom_cmap,#cmap='YlOrRd',#cmap=custom_cmap, # 使用自定义的颜色映射legend=False,#不使用默认的比例尺legend_kwds={'label': "Insurance coverage rate(%)",'orientation': "horizontal", # 设置比例尺为水平'pad': 0.005, # 调整比例尺与图形之间的间距},missing_kwds={'color': 'none', # 设置无数据区域的背景为透明'hatch': '///', # 无数据区域使用斜黑线填充'label': 'No Data' # 在图例中添加说明},edgecolor='black',ax=ax
)# 去掉地图的经纬度尺度
ax.axis('off')# 获取 colorbar 对象并调整其位置
cbar = fig.colorbar(ax.collections[0], ax=ax, orientation='horizontal') # 获取colorbar对象
cbar.ax.set_position([0.25, 0.05, 0.5, 0.03]) # 设置colorbar的位置# 获取 colorbar 的 Axes 对象，调整比例尺数字和标签
cbar.ax.xaxis.set_ticks_position('top') # 设置刻度到上方
cbar.ax.xaxis.set_label_position('top') # 设置标签到上方
cbar.ax.set_xlabel("Insurance coverage rate(%)", labelpad=10) # 设置标签，并调整与比例尺的间距# 添加无数据方块到比例尺左边
no_data_rect = Rectangle((0.1, 0.05), 0.05, 0.03, transform=fig.transFigure, facecolor='none', edgecolor='black', hatch='///')
fig.patches.append(no_data_rect)# 添加 "No Data" 标注文字到比例尺上方
ax.text(0.125, 0.1, "No Data", transform=fig.transFigure, fontsize=10, verticalalignment='center', horizontalalignment='center')plt.show()
地图js文件
文件过大，可以去GeoJSON Maps of the globe等诸多网页上下载，或者找我我私发
滞后图
import numpy as np
import matplotlib.pyplot as plt# 上层食物限制系数
aP = 0.002
aI = 0.004
aI2 = 0.004
aS = 0.005
aBat = 0.005
aBird = 0.005# 捕食系数和a一个数量级
bI = 0.002
bI2 = 0.002
bS = 0.001
bBat = 0.002
bBird = 0.001
bF = 0.001# 竞争系数
cP = 0.018
cG = 0.02
cI = 0.007
cI2 = 0.005
cBat = 0.004
cBird = 0.004# 自然死亡
dI, dI2, dS = 0.2, 0.18, 0.07
dBat, dBird = 0.08, 0.06
dF = 0.05# 化学影响
wa = 1
wb = 0.5# 初始数量
P0, G0, I0, I20, S0 = 10, 5, 2, 3, 1
Bat0, Bird0 = 1, 1
F0 = 1# 自然增长
rP, rG, rI, rI2, rS = 3, 3.3, 2.9, 2.8, 2.5
rBat, rBird = 2.6, 2.4
rF = 2# 限制
kP, kG, kI, kI2, kS = 100, 50, 20, 25, 10
kBat, kBird = 10, 8
kF = 5# 设置时长和最大时间
dt = 0.1
T = 12
time_steps = int(T / dt)# 初始化数量数组
P = np.zeros(time_steps)
G = np.zeros(time_steps)
I = np.zeros(time_steps)
I2 = np.zeros(time_steps)
S = np.zeros(time_steps)
F = np.zeros(time_steps)# 初始值
P[0], G[0], I[0], I2[0], S[0], F[0] = P0, G0, I0, I20, S0, F0# 动态更新公式
def P0_t(t):return rP * (1 - P[t-1] / kP) * P[t-1] - bI * I[t-1] * P[t-1]def G0_t(t):return rG * (1 - G[t-1] / kG) * G[t-1]def I0_t(t):return rI * (1 - I[t-1] / kI) * I[t-1] - dI * I[t-1] - bS * S[t-1] * I[t-1] - aP * P[t-1] * I[t-1]def S0_t(t):return rS * (1 - S[t-1] / kS) * S[t-1] - dS * S[t-1] - aI * I[t-1] * S[t-1]# 动态模拟
for t in range(1, time_steps):DP = P0_t(t) - bI2 * I2[t-1] * P[t-1]DG = G0_t(t) - cP * P[t-1] * G[t-1] - wa * G[t-1]DI = I0_t(t) - wb * I[t-1] - bBat * I[t-1]DI2 = rI2 * (1 - I2[t-1] / kI2) * I2[t-1] - dI2 * I2[t-1]DS = S0_t(t) - aI * I[t-1] * S[t-1]DF = rF * (1 - F[t-1] / kF) - dF * F[t-1]P[t] = P[t-1] + DP * dt / 5G[t] = G[t-1] + DG * dt / 5I[t] = I[t-1] + DI * dt / 5I2[t] = I2[t-1] + DI2 * dt / 5S[t] = S[t-1] + DS * dt / 5F[t] = F[t-1] + DF * dt / 5# 滞后图绘制
variables = {"Crops": P,"Weeds": G,"Insects": I,"Voles": I2,"Secondary consumers": S,"Tertiary consumers": F
}
'''
# 自定义颜色字典
custom_colors = {"Crops (P)": "#D2A000","Weeds (G)": "#8E6E3A","Insects (I)": "#647B35","Voles (I2)": "#56766E","Secondary consumers (S)": "#2C7488","Tertiary consumers (F)": "#85352F"
}
'''
custom_colors = {"Crops": "#FFCA21","Weeds": "#D5BE98","Insects": "#A2BE6A","Voles": "#8BABA3","Secondary consumers": "#62B5CC","Tertiary consumers": "#E5B9B5"
}fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()for i, (label, data) in enumerate(variables.items()):ax = axes[i]ax.plot(data[:-1], data[1:], marker='o', linewidth=2, label=label, color=custom_colors[label])ax.scatter(data[0], data[1], color="green", label="Start", s=80)ax.scatter(data[-2], data[-1], color="red", label="End", s=80)ax.set_xlabel('$y(t)$')ax.set_ylabel('$y(t+1)$')ax.set_title(f'Lag Diagram: {label}')ax.legend()ax.grid(True)plt.tight_layout()
plt.show()
斯皮尔曼相关矩阵
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap# 直接输入相关性矩阵
corr_matrix = np.array([[1, 5/7, 5/3, 1, 5, 5/3, 5, 5/3],[7/5, 1, 7/3, 7/5, 7, 7/3, 7, 7/3],[3/5, 3/7, 1, 3/5, 3, 1, 3, 1],[1, 5/7, 5/3, 1, 5, 5/3, 5, 5/3],[1/5, 1/7, 3, 1/5, 1, 1/3, 1, 1/3],[3/5, 3/7, 1, 3/5, 3, 1, 3, 1],[1/5, 1/7, 1/3, 1/5, 1, 1/3, 1, 1/3],[3/5, 3/7, 1, 3/5, 3, 1, 3, 1]
])# 转换为 DataFrame，提供正确的行列名
corr_matrix = pd.DataFrame(corr_matrix,columns=['1','2','3','4','5','6','7','8'], # 8 个变量列名index=['1 - PH','2 - Organic Substance','3 - Water Pollution','4 - Chemical Residues','5 - Biological Tiers','6 - Diversity','7 - Cost','8 - Revenue'] # 8 个变量行名
)# 创建掩码矩阵
mask_upper = np.triu(np.ones_like(corr_matrix, dtype=bool), k=0) # 上三角掩码（包括对角线）
mask_lower = ~mask_upper # 下三角排掩码（不包括对角线）# 设置左下角绿色背景矩阵
left_lower_colors = corr_matrix.applymap(lambda x: "#FDD57F" if x >= 3 or x <= 1/3 else "white"
)# 设置绘图风格
plt.figure(figsize=(10, 8))
sns.set_theme(style="white")# 自定义配色：红色到蓝色的渐变
custom_cmap = LinearSegmentedColormap.from_list("red_blue", ["red", "white", "#8596C4"], N=100 # 红色到蓝色渐变，设置更多的渐变色
)# 绘制右上三角（颜色渐变，带标记）
sns.heatmap(corr_matrix,mask=mask_lower, # 掩盖左下部分cmap=custom_cmap, # 使用自定义颜色annot=False, # 不显示原始数据值vmin=0, vmax=7, # 设置颜色范围cbar_kws={"shrink": 0.8}, # 调整颜色条大小linewidths=0.5, # 方块分割线宽度linecolor="black", # 分割线颜色为黑色square=True, # 方形格子cbar=True # 显示颜色条
)# 绘制左下三角（填充背景色，显示数值，并保留分隔线）
for i in range(corr_matrix.shape[0]):for j in range(i): # 左下三角plt.gca().add_patch(plt.Rectangle((j, i), 1, 1, facecolor=left_lower_colors.iloc[i, j], edgecolor="black", linewidth=0.5))plt.text(j + 0.5, # X 坐标i + 0.5, # Y 坐标f"{corr_matrix.iloc[i, j]:.2f}", # 显示数值（两位小数）ha="center", va="center", # 水平和垂直居中color="black", # 数值字体颜色fontsize=14, # 数值字体大小)# 绘制右上三角的圆圈符号（满足条件时显示）
for i in range(corr_matrix.shape[0]):for j in range(i + 1, corr_matrix.shape[1]): # 右上三角if corr_matrix.iloc[i, j] >= 3 or corr_matrix.iloc[i, j] <= 1/3: # 满足条件plt.text(j + 0.5, # X 坐标i + 0.5, # Y 坐标"o", # 使用圆圈符号ha="center", va="center", # 水平和垂直居中color="#9D0109", # 圆圈符号颜色fontsize=18, # 符号字体大小fontweight="bold" # 加粗字体)# 添加标题
plt.title("Spearman Matrix Heatmap", fontsize=16)
plt.show()
拟合曲线
import numpy as np
import matplotlib.pyplot as plt# 定义 x 轴范围
x = np.linspace(-7, 50, 500)# 分段点
x0, x1, x2, x3 = -1, 3, 7, 12def curve0(x):return 3.25 + 4 * np.log(-x)# 定义每段曲线
def curve1(x):y_start = curve0(x0)y_end = -0.2# 第一段：下凹地增加return -2 * (1 - np.exp(-x))def curve2(x):# 第二段：下凸地减小y_start = curve1(x1)y_end = -0.4return y_start - (y_start - y_end) * ((x - x1) / (x2 - x1))**2def curve3(x):# 第三段：下凸地增加y_start = curve2(x2)y_end = -0.6return -8.2 + 4 * np.log(x)def curve4(x):return 5.9 - 50 / x# 拼接曲线
y = np.piecewise(x,[x < x0, (x >= x0) & (x < x1), (x >= x1) & (x < x2), (x > x2) & (x <= x3), x > x3],[curve0, curve1, curve2, curve3, curve4]
)# 绘图
plt.figure(figsize=(10, 6))
plt.plot(x, y, label='kkk', color='#578B38',linewidth=5)# 去除x轴的数字
plt.xticks([])
plt.yticks([])# 去除网格线
plt.grid(False)
# 去除边框
ax = plt.gca()
ax.spines['top'].set_visible(False) # 去除顶部边框
ax.spines['right'].set_visible(False) # 去除右侧边框
ax.spines['left'].set_visible(False) # 去除左侧边框
ax.spines['bottom'].set_visible(False) # 去除底部边框
'''
# 设置边框（保留边框）
ax = plt.gca()
ax.spines['bottom'].set_linewidth(1) # 设置x轴轴线宽度
ax.spines['left'].set_linewidth(1) # 设置y轴轴线宽度
'''
'''
# 自定义y轴刻度
plt.yticks([300, 500, 700, 900]) # 自定义y轴刻度值
'''
'''
# 设置背景颜色
ax.axvspan(-7, x0, color='#C4D6A0', alpha=0.3)
ax.axvspan(x0, x2, color='#F6F8EB', alpha=0.3) # 在 x0 到 x1 范围内填充黄色背景，透明度为 0.3
ax.axvspan(x2, 25, color='#FFFF9E', alpha=0.5) # 在 x2 到 x3 范围内填充浅蓝色背景，透明度为 0.5
#ax.axvspan(x1, x2, color='lightgreen', alpha=0.3) # 在 x0 到 x1 范围内填充黄色背景，透明度为 0.3
ax.axvspan(25, 50, color='#D5BE98', alpha=0.5) # 在 x2 到 x3 范围内填充浅蓝色背景，透明度为 0.5
'''
# 使用 fill_between 填充不同区间
plt.fill_between(x[x < x0], y[x < x0], -5, color='#FFFF9E', alpha=0.3, label='Region 1')
plt.fill_between(x[(x >= x0) & (x < x2)], y[(x >= x0) & (x < x2)], -5, color='#C7EC99', alpha=0.3, label='Region 2')
plt.fill_between(x[(x >= x2) & (x <= 25)], y[(x >= x2) & (x <= 25)], -5, color='#8ED694', alpha=0.5, label='Region 3')
plt.fill_between(x[x > 25], y[x > 25], -5, color='#36B88E', alpha=0.5, label='Region 4')# 去除左右两边的留白
plt.xlim(np.min(x), np.max(x))# 显示图形
plt.show()
敏感性分析
import numpy as np
import matplotlib.pyplot as plt# 定义公式参数
N0 = 10 # 初始值
r = 0.1 # 增长率
K = 200 # 环境容量
T = 3/2 # 周期
time = np.linspace(0, 3.5, 1000) # 时间范围
epsilon = np.random.normal(0, 0.5, len(time)) # 噪声# 定义函数 N(t)
def N(t, N0, r, K, T, epsilon):result = N0 * np.exp(r * (1 - N0 / K)) * np.sin(2 * np.pi * t / T) + epsilon'''if t<1:result = N0 * np.exp(r * (1 - N0 / K)) * np.sin(2 * np.pi * t / T) + epsilonelif t>1.5:return'''result[result <= 0] = np.nan # 过滤掉 y <= 0 的部分return result*15# 计算 N(t)
N_values = N(time, N0, r, K, T, epsilon)# 绘制敏感性分析（调整参数 r 和 K）
plt.figure(figsize=(10, 6))# r 敏感性分析
for r_test in [0.05, 0.1, 0.2]:N_values_r = N(time, N0, r_test, K, T, epsilon)plt.plot(time, N_values_r, label=f"r={r_test*30}")# K 敏感性分析
for K_test in [100, 200, 400]:N_values_K = N(time, N0, r, K_test, T, epsilon)plt.plot(time, N_values_K, linestyle="--", label=f"K={K_test}")plt.title("Sensitivity Analysis of Agricultural Seasonality")
plt.xlabel("Time")
plt.ylabel("Number of Crops")
plt.legend()
plt.grid()
plt.show()