如何在地图（底图）上绘制插值站数据？

Question

我有一个文件，其中包含 5 个站点超过一个月的累积降雨量。 csv 文件中有纬度、经度和降雨数据。我的文件是这样的：

Out[18]: 
         lat       lon  rain
0 -48.379000 -1.067000  213.0
1 -48.435548 -1.401513  157.2
2 -48.482217 -1.449707  147.0
3 -48.457779 -1.249272  182.6
4 -48.479847 -1.308735   49.4

我正在努力：

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

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)

 data = pd.read_csv('.../rainfall-2010-12.txt',na_values=['NaN'], sep=',')
norm = Normalize()

#mapextent
lllon = data['lon'].min()
lllat = data['lat'].min()
urlon = data['lon'].max()
urlat = data['lat'].max()

#Basemap
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))

#griddata
numcols, numrows = 300, 300
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.rain.values
zi = griddata(x, y, z, xi, yi, interp='linear')

m.drawcoastlines()

# contour plot
conf = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')

cbar = plt.colorbar(conf, orientation='horizontal', fraction=.057, pad=0.05)
cbar.set_label("Rainfall - mm")

plt.title("Rainfall")
plt.show()

但是当我尝试运行时，我收到了这个错误消息：

/usr/local/lib/python3.5/dist-packages/mpl_toolkits/basemap/__init__.py:3608: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.   b = ax.ishold() /usr/local/lib/python3.5/dist-packages/mpl_toolkits/basemap/__init__.py:3675: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.   ax.hold(b) Traceback (most recent call last):

  File "<ipython-input-17-cb8133160e02>", line 4, in <module>
    conf = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')

  File "/usr/local/lib/python3.5/dist-packages/mpl_toolkits/basemap/__init__.py", line 521, in with_transform
    return plotfunc(self,x,y,data,*args,**kwargs)

  File "/usr/local/lib/python3.5/dist-packages/mpl_toolkits/basemap/__init__.py", line 3644, in contourf
    xx = x[x.shape[0]/2,:]

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

我该如何解决这个问题？

Answer 1

您的代码中没有严重错误。 (1) 我怀疑 (Lat, long) 是颠倒的，因为所有点都位于南大西洋。 (2) 数据文件内容第一行不能有多余的空格。

这是数据。第一行没有多余的空格。

lat,lon,rain
0, -48.379000, -1.067000,  213.0
1, -48.435548, -1.401513,  157.2
2, -48.482217, -1.449707,  147.0
3, -48.457779, -1.249272,  182.6
4, -48.479847, -1.308735,   49.4

这是基于您的工作代码和生成的地图。 请注意，我注释掉了绘制土地特征的部分。它们对当前数据毫无用处。

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 math import ceil

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)

data = pd.read_csv('rainfall-2010-12.txt', na_values=['NaN'], sep=',')
norm = Normalize()

#mapextent
lllon = data['lon'].min()
lllat = data['lat'].min()
urlon = data['lon'].max()
urlat = data['lat'].max()

#Basemap
pad = 0.01  # padding around map extents
m = Basemap(
    projection = 'merc', \
    llcrnrlon = lllon - pad, \
    llcrnrlat = lllat - pad, \
    urcrnrlon = urlon + pad, \
    urcrnrlat = urlat + pad, \
    resolution='h', \
    ax=ax)

# transform lon / lat coordinates to map projection
data['projected_lon'], data['projected_lat'] = m(*(data.lon.values, data.lat.values))

#griddata
numcols, numrows = 100, 100
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.rain.values
zi = griddata(x, y, z, xi, yi, interp='linear')

# contours plot
conf = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu')  # filled contour
cont = m.contour(xi, yi, zi, zorder=5)                           # line contour

# *** If the location is on land, uncomment this block of code ***
# draw some map features
#m.drawcoastlines()
#m.drawrivers(color='b')
#m.fillcontinents(color='lightyellow')

# paint the ocean
watercolor=np.array([0.6, 0.8, 0.9])
m.drawlsmask(land_color=watercolor, ocean_color=watercolor, lakes=False, resolution='i')

# draw parallels and meridians; also accompanying labels
m.drawparallels(np.arange(ceil(lllat*100)/100, ceil(urlat*100)/100, 0.05), \
                labels=[1,0,0,0], linewidth=0.2, dashes=[5, 3])   # draw parallels
m.drawmeridians(np.arange(ceil(lllon*100)/100, ceil(urlon*100)/100, 0.05), \
                labels=[0,0,0,1], linewidth=0.2, dashes=[5, 3])   # draw meridians

cbar = plt.colorbar(conf, orientation='horizontal', fraction=.057, pad=0.075)
cbar.set_label("Rainfall - mm")

plt.title("Rainfall")
plt.show()

（编辑 1） 使用上面相同的代码，但有两种不同的设置：

(1) Python 2.7.14 / 底图 1.1.0 (conda-forge)

(2) Python 3.5.4 / 底图 1.1.0 (conda-forge)

生成的图在视觉上与下图相同， 左图：设置 1，右图：设置 2。