利用Python插值繪製等值線圖
阿新 • • 發佈:2018-12-04
最近需要根據有限的站位點繪製插值等值線圖,在網上中文搜尋一通,只發現了這貨Matplot Basemap 畫湖北地圖、插值、等值線,要麼就是對這貨的轉載,這貨不提供資料的形式,但是基本的程式碼思路還是不錯的,於是繼續轉向google英文,搜到瞭如下的回答,我對程式碼做了註釋,已備別人查詢,關於文中提到的資料是txt格式的,我也直接將資料貼在下面了。總結一下:在地圖上繪製等值線:
- 確定基本的繪圖框架;
- 獲取採集資料,與地圖做對映,並根據對映資料插值;scipy.interpolate.griddata包插值比較快,常用的三種插值方法為liner(基於三角形的線性插補法),cubic(基於三角形的三次插補法),nearest( 最近鄰居插補法),這些方法定義了匹配資料點的曲面型別,'cubic' 方法生成平滑曲面,而 'linear' 和 'nearest' 分別具有一階導數和零階導數不連續。
- 根據柵格插值資料繪圖
# -*- coding: utf-8 -*- @author: Adwiy Wang import numpy as np import pandas as pd from matplotlib.mlab import griddata from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt from matplotlib.colors import Normalize from scipy.interpolate import griddata as gd # 設定基本圖片畫板 fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, axisbg='w', frame_on=False) # 提取資料 data = pd.read_csv('datam.txt', delim_whitespace=True) norm = Normalize() #設定地圖邊界值 lllon = 21 lllat = -18 urlon = 34 urlat = -8 #初始化地圖 m = Basemap( projection = 'merc', llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat, resolution='h') # 將經緯度點轉換為地圖對映點 data['projected_lon'], data['projected_lat'] = m(*(data.Lon.values, data.Lat.values)) # 生成經緯度的柵格資料 numcols, numrows = 1000, 1000 xi = np.linspace(data['projected_lon'].min(), data['projected_lon'].max(), numcols) yi = np.linspace(data['projected_lat'].min(), data['projected_lat'].max(), numrows) xi, yi = np.meshgrid(xi, yi) # 插值 x, y, z = data['projected_lon'].values, data['projected_lat'].values, data.Z.values zi = gd( (data[['projected_lon', 'projected_lat']]), data.Z.values, (xi, yi), method='cubic') # 設定地圖細節 m.drawmapboundary(fill_color = 'white') m.fillcontinents(color='#C0C0C0', lake_color='#7093DB') m.drawcountries( linewidth=.75, linestyle='solid', color='#000073', antialiased=True, ax=ax, zorder=3) m.drawparallels( np.arange(lllat, urlat, 2.), color = 'black', linewidth = 0.5, labels=[True, False, False, False]) m.drawmeridians( np.arange(lllon, urlon, 2.), color = '0.25', linewidth = 0.5, labels=[False, False, False, True]) # 等值面圖繪製 con = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='jet') # 插入測繪點 m.scatter( data['projected_lon'], data['projected_lat'], color='#545454', edgecolor='#ffffff', alpha=.75, s=50 * norm(data['Z']), cmap='jet', ax=ax, vmin=zi.min(), vmax=zi.max(), zorder=4) # 插入色標、名稱和範圍 cbar = plt.colorbar(con,orientation='horizontal', fraction=.057, pad=0.05) cbar.set_label("Mean Rainfall - mm") m.drawmapscale( 24., -9., 28., -13, 100, units='km', fontsize=10, yoffset=None, barstyle='fancy', labelstyle='simple', fillcolor1='w', fillcolor2='#000000', fontcolor='#000000', zorder=5) plt.title("Mean Rainfall") plt.savefig("rainfall.png", format="png", dpi=300, transparent=True) plt.show()
資料檔案:datam.txt
on Lat Z 32.6 -13.6 41 27.1 -16.9 43 32.7 -10.2 46 24.2 -13.6 33 28.5 -14.4 43 28.1 -12.6 33 27.9 -15.8 46 24.8 -14.8 44 31.1 -10.2 35 25.9 -13.5 24 29.1 -9.8 10 25.8 -17.8 39 33.2 -12.3 44 28.3 -15.4 46 27.6 -16.1 47 28.9 -11.1 31 31.3 -8.9 39 31.9 -13.3 45 23.1 -15.3 31 31.4 -11.9 39 27.1 -15.0 42 24.4 -11.8 15 28.6 -13.0 39 31.3 -14.3 44 23.3 -16.1 39 30.2 -13.2 38 24.3 -17.5 32 26.4 -12.2 23 23.1 -13.5 27