我几个月前问过这个问题,但仍在努力寻求解决方案。我的代码并排提供了底图和轮廓图(但是打印到文件仅给出轮廓图),但是我希望它们叠加在一起。最好的解决方案将是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
You're almost there, but Basemap can be temperamental, and you have to manage the z-order of plots / map details. Also, you have to transform your lon / lat coordinates to map projection coordinates before you plot them using basemap.
Here's a complete solution, which gives the following output. I've changed some colours and linewidths in order to make the whole thing more legible, YMMV. I've also scaled the size of the scatter points by the normalised 'mean' value (data['Z']
) – you can simply remove it and substitute e.g. 50
if you'd prefer a constant size (it'll look like the largest marker).
Please also specify the units of rainfall, and the duration of the measurement which generated the means, if possible:
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
)提供最好的产生阴谋。
我还要补充一点,上面演示和讨论的插值方法是最简单的方法,不一定对降雨数据的插值有效。本文对用于水文学和水文建模的有效方法和注意事项进行了很好的概述。这些的实现(可能使用Scipy)留作练习&c。
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句