# 準(zhǔn)備數(shù)據(jù) 讀取數(shù)據(jù)
# 0.讀寫(xiě)實(shí)際數(shù)據(jù)生成三維曲面，數(shù)據(jù)格式為x y z
filename=r'D:\project\PythonProject\ECL\data\geochemical-data\2018_T28.txt'
dataTop = pd.read_csv(filename, sep='\t', header=None, names=['x', 'y', 'z'])
# 去掉無(wú)效數(shù)據(jù)，一般為-99999.0000
data = dataTop[dataTop['z'] != -99999.0000]
x = data.iloc[:, 0]
y = data.iloc[:, 1]
z = data.iloc[:, 2]*(-1) #深度值是負(fù)數(shù)，要取反。
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
xi, yi = np.meshgrid(xi, yi) # 將一維數(shù)據(jù)處理為二維的網(wǎng)格數(shù)據(jù)
zi = griddata(data.iloc[:,0:2], z, (xi, yi), method='cubic') # 用法詳見(jiàn)附錄1

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(x, y, z, color='white', edgecolors='grey', alpha=0.5) #繪制三角網(wǎng)格組成的三維曲面
ax.scatter(x, y, z, c='red')  # 繪制三維散點(diǎn)圖
plt.show()

fig = plt.figure()
# ax = plt.axes(projection='3d')
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xi, yi, zi, cmap='BuPu', linewidth=0, antialiased=True)  #繪制三維曲面
# surf = ax.scatter(xi, yi, zi, cmap='BuPu', linewidth=0, antialiased=True)     #繪制三維散點(diǎn)圖
# surf = ax.contourf(xi, yi, zi, zdim='z',offset=0.3, cmap='BuPu')   #等高線面圖（contourf）或等高線圖（contour），要設(shè)置offset，為Z的最小值，
fig.colorbar(surf)
ax.set_title('三維圖')
ax.set_zlim3d(np.min(z), np.max(z))
plt.show()

filename = r'D:\project\PythonProject\ECL\data\geochemical-data\0.txt'  # 數(shù)據(jù)文件地址，附件1
df = pd.read_csv(filename, sep="\t")  # 讀取文件
df1 = df["1"]  # 讀取第一列數(shù)據(jù)
df2 = df['2']  # 讀取第二列數(shù)據(jù)
df3 = df['3']  # 讀取第三列數(shù)據(jù)
odf1 = np.linspace(100, 1900, 50)  # 設(shè)置網(wǎng)格經(jīng)度
odf2 = np.linspace(10, 600, 50)  # 設(shè)置網(wǎng)格緯度
odf1, odf2 = np.meshgrid(odf1, odf2)  # 網(wǎng)格化，生成網(wǎng)絡(luò),生成網(wǎng)格形狀是第一個(gè)維度對(duì)應(yīng)odf1，第二個(gè)維度對(duì)應(yīng)odf2
func = Rbf(df1, df2,df3, function='linear')  # 定義插值函數(shù)plt.cm.hot
odf3_new = func(odf1, odf2)  # 獲得插值后的網(wǎng)格累計(jì)降水量
fig = plt.figure(figsize=(12, 7))
ax1 = plt.axes(projection='3d')  # 創(chuàng)建三維坐標(biāo)軸
ax1.plot_surface(odf1, odf2, odf3_new,alpha=0.3,cmap='rainbow') #繪制三維曲面,alpha-控制透明度,cmap-控制顏色
# 繪制z方向投影填充圖，等高線面圖，投到x-y平面，offset為z最小值。
cs=plt.contourf(odf1, odf2, odf3_new,zdir='z',offset=0,
             levels=np.arange(odf3_new.min(), odf3_new.max(), (odf3_new.max() - odf3_new.min()) / 10), cmap='GnBu',
             extend='both')  # 畫(huà)圖
# 繪制等高線圖
line = plt.contour(odf1, odf2, odf3_new,zdir='z',offset=0,cmap="rainbow",levels=np.arange(odf3_new.min(), odf3_new.max(), (odf3_new.max() - odf3_new.min()) / 10))
plt.clabel(line, inline=True, fontsize=12)
ax1.set_title('降雨量三維等值線圖')
plt.colorbar(cs)
plt.show()

# 繪制二維等值線
levels = np.linspace(np.min(z), np.max(z), 50)
fig, ax = plt.subplots(figsize=(8, 6))
# # 1.設(shè)置顏色條  第一種方式
# print(cm.colors.cnames)
# 默認(rèn)顏色標(biāo)如 jet，coolwarm，gnuplot2_r,RdBu_r,PuBuGn_r,ocean_r，輸入的顏色名稱(chēng)錯(cuò)誤時(shí)，會(huì)自動(dòng)輸出色標(biāo)的列表
cmap = cm.get_cmap('seismic_r')
# cmap = cm.get_cmap('jet', 10)  # 將色條分成10截
norm = cm.colors.Normalize(vmin=np.min(z), vmax=np.max(z))  # 設(shè)置色條表示的數(shù)值范圍
im1 = cm.ScalarMappable(norm=norm, cmap=cmap)               # 設(shè)置映射很重要
# # 繪制顏色條(left, bottom, width, height)--表示figure的百分比,從figure 從橫向92%，縱向10%的位置開(kāi)始繪制, 寬是figure的3%，高是figure的78%,
ax9 = fig.add_axes([0.92, 0.1, 0.03, 0.78])
cb = plt.colorbar(im1, cax=ax9, orientation='vertical', extend='neither') #縱向繪制，兩端無(wú)箭頭
# ticks與norm對(duì)應(yīng)
# # cb = plt.colorbar(im1, cax=ax9, orientation='horizontal', extend='max', ticks=np.linspace(1900,2600, 51))
# cs = ax.contour(xi, yi, zi, levels=levels, cmap=cmap)         # 不存在顏色間隔分段，并指定顏色條
# cs = ax.contour(xi, yi, zi, levels=levels,cmap='coolwarm')  # 存在顏色間隔分段
cs = ax.contourf(xi, yi, zi, levels=levels,cmap='jet',extend='neither') # 等值線填充
ax.clabel(cs, inline=True, fontsize=6)
ax.set_title('等高線圖')
plt.show()
# # 2.設(shè)置顏色條 第二種方式
# cs = ax.contourf(xi, yi, zi, levels=levels, cmap='jet', extend='neither')  # 等值線填充，存在顏色間隔分段
# # ax.clabel(cs, inline=True, fontsize=6)
# ax.set_title('等高線圖')
# plt.colorbar(cs)
# plt.show()

# 1.矢量場(chǎng)流線圖 樣例
# 0:5表示數(shù)組中數(shù)值所在的區(qū)間。100j表示劃分的密度，值越大，圖片越清晰
Y1, X1 = np.mgrid[-5:5:1000j, -5:5:10j]
# （X,Y）是一維numpy數(shù)組的等距網(wǎng)格，（U,V）參數(shù)匹配的是（X,Y）速率的二維numpy數(shù)組
# U,V矩陣在維度上的行數(shù)必須等于Y的長(zhǎng)度，列的數(shù)量必須匹配X的長(zhǎng)度
U = -1 - X1 ** 2 + Y1
V = 1 + X1 - Y1 ** 2
# 可視化矢量場(chǎng)
# 矢量場(chǎng)中的種子點(diǎn)坐標(biāo)
seed_points = np.array([[-2, -1, 0, 1, 2, -1], [-2, -1, 0, 1, 2, 2]]) # 種子點(diǎn)
# cs=plt.streamplot(X1, Y1, U, V,density=[0.5,1],color=U,cmap="autumn",linewidth=1,start_points=seed_points.T)
cs = plt.streamplot(X1, Y1, U, V, color=U, cmap="Accent", linewidth=1)
# plt.plot(seed_points[0], seed_points[1], "+", color="g") # 繪制折線圖，使用marker屬性標(biāo)記
plt.colorbar()
plt.show()
# 2.用實(shí)際數(shù)據(jù)進(jìn)行繪制,UV如何計(jì)算得到，要根據(jù)不同的目的進(jìn)行計(jì)算。
U,V=vectorComputeUV(xi,yi,zi)
cs = plt.streamplot(xi, yi, U, V, color=U, cmap="Accent", linewidth=1)
# cs = plt.quiver(xi, yi, U, V)
plt.colorbar()
plt.show()
# 計(jì)算矢量場(chǎng)的速度矢量
def vectorComputeUV(xi,yi,zi):
    # U = np.log10((xi/10000)) - zi/1000
    # V = 1 + np.log10(yi/100000) - (zi/1000) ** 3
    U = -1 - (xi/10000) ** 2 + (yi/100000)
    V = 1 + (xi/10000)  - (yi/100000) ** 2
    return U,V