这是我几个月前提出的问题,现在仍在努力寻找解决方案。我的代码并排给我一个 basemap 和一个等高线图(但打印到文件只给出等高线图),但我希望它们叠加。最好的解决方案是这里的 https://gist.github.com/oblakeobjet/7546272但这并没有说明如何引入数据,在线从头学习时很难。 在不累非常善良的人的情况下,我希望解决方案就像更改一行代码一样简单,并且有人可以提供帮助。 我的代码
#!/usr/bin/python
# vim: set fileencoding=UTF8
import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
#fig = plt.figure(figsize=(10,8)) #when uncommented draws map with colorbar but no contours
#prepare a basemap
m = Basemap(projection = 'merc',llcrnrlon = 21, llcrnrlat = -18, urcrnrlon = 34, urcrnrlat = -8, resolution='h')
# draw country outlines.
m.drawcountries(linewidth=0.5, linestyle='solid', color='k', antialiased=1, ax=None, zorder=None)
m.drawmapboundary(fill_color = 'white')
m.fillcontinents(color='coral',lake_color='blue')
parallels = np.arange(-18, -8, 2.)
m.drawparallels(np.arange(-18, -8, 2.), color = 'black', linewidth = 0.5)
m.drawparallels(parallels,labels=[True,False,False,False])
meridians = np.arange(22,34, 2.)
m.drawmeridians(np.arange(21,36, 2.), color = '0.25', linewidth = 0.5)
m.drawmeridians(meridians,labels=[False,False,False,True])
fig = plt.figure(figsize=(10,8)) # At this position or commented draws teo figures side by side
#-- Read the data.
data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True)
#-- Now gridding data. First making a regular grid to interpolate onto
numcols, numrows = 300, 300
xi = np.linspace(data.Lon.min(), data.Lon.max(), numcols)
yi = np.linspace(data.Lat.min(), data.Lat.max(), numrows)
xi, yi = np.meshgrid(xi, yi)
#-- Interpolating at the points in xi, yi
x, y, z = data.Lon.values, data.Lat.values, data.Z.values
zi = griddata(x, y, z, xi, yi)
#-- Display and write the results
m = plt.contourf(xi, yi, zi)
plt.scatter(data.Lon, data.Lat, c=data.Z, s=100,
vmin=zi.min(), vmax=zi.max())
fig.colorbar(m)
plt.savefig("rainfall.jpg", format="jpg")
我得到的图看起来像这样 和
和我的数据
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
最佳答案
您就快完成了,但是 basemap 可能不稳定,您必须管理绘图/ map 细节的 z 顺序。此外,在使用 basemap 绘制它们之前,您必须将经/纬度坐标转换为 map 投影坐标。
这是一个完整的解决方案,它提供以下输出。我更改了一些颜色和线宽,以使整个内容更清晰易读,YMMV。我还通过标准化的“平均”值 (data['Z']
) 缩放了散点的大小——您可以简单地删除它并替换,例如50
如果您更喜欢恒定大小(它看起来像最大的标记)。
如果可能的话,还请具体说明降雨量的单位和生成平均值的测量持续时间:
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
%matplotlib inline
# set up plot
plt.clf()
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, axisbg='w', frame_on=False)
# grab data
data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True)
norm = Normalize()
# define map extent
lllon = 21
lllat = -18
urlon = 34
urlat = -8
# Set up Basemap instance
m = Basemap(
projection = 'merc',
llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat,
resolution='h')
# transform lon / lat coordinates to map projection
data['projected_lon'], data['projected_lat'] = m(*(data.Lon.values, data.Lat.values))
# grid data
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)
# interpolate
x, y, z = data['projected_lon'].values, data['projected_lat'].values, data.Z.values
zi = griddata(x, y, z, xi, yi)
# draw map details
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])
# contour plot
con = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')
# scatter plot
m.scatter(
data['projected_lon'],
data['projected_lat'],
color='#545454',
edgecolor='#ffffff',
alpha=.75,
s=50 * norm(data['Z']),
cmap='RdPu',
ax=ax,
vmin=zi.min(), vmax=zi.max(), zorder=4)
# add colour bar and title
# add colour bar, title, and scale
cbar = plt.colorbar(conf, 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()
使用 matplotlib 的 griddata
方法很方便,但也可能很慢。作为替代方案,您可以使用 scipy 的 griddata更快、更灵活的方法:
from scipy.interpolate import griddata as gd
zi = gd(
(data[['projected_lon', 'projected_lat']]),
data.Z.values,
(xi, yi),
method='linear')
如果您使用 scipy 的 griddata
方法,您还必须确定哪个方法(nearest
、linear
、 cubic
) 给出最好的结果图。
我应该补充一点,上面演示和讨论的插值方法是最简单的方法,不一定对降雨数据的插值有效。 This article很好地概述了在水文学和水文建模中使用的有效方法和注意事项。这些的实现(可能使用 Scipy)留作练习 &c。
关于python - 如何将等高线图叠加在 basemap 上,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/26872337/