如何在地图(底图)上绘制插值站数据? [英] How to plot interpolate station data over a map (basemap)?

问题描述

我有一个文件，该文件在5个站点中累积了一个月的降雨量. csv文件中有纬度，经度和降雨数据.我的文件就是这样:

I've got a file with accumulated rainfall over a month in 5 stations. There are lat, lon and rain data in csv file. My file is just like that:

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()

但是当我尝试运行时，出现了此错误消息:

But when I'm trying to run, I got this error msg:

/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

我该如何解决?

推荐答案

您的代码中没有严重错误. (1)由于所有点都位于南大西洋，我怀疑(纬长)是相反的. (2)数据文件中的内容第一行不应有多余的空格.

There is no serious error in your code. (1) I suspect that (Lat, long) is reversed because all points are located in South Atlantic. (2) Content in data file should not have extra spaces on the first line.

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

This is the data. No extra spaces on the first line.

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

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

Here is the working code based on yours, and the resulting map. Note that I commented out the parts that plot land features. They are useless with the current data.

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) 上面的代码相同，但是有2种不同的设置:

(Edit 1) With the same code above, but 2 different setups:

(1)Python 2.7.14/Basemap 1.1.0(conda-forge)

(1) Python 2.7.14 / Basemap 1.1.0 (conda-forge)

(2)Python 3.5.4/Basemap 1.1.0(conda-forge)

(2) Python 3.5.4 / Basemap 1.1.0 (conda-forge)

结果图在视觉上与以下所示相同， 左图:设置1，右图:设置2.

The resulting plots are visually the same as shown below, Left image: setup 1, right image: setup 2.

这篇关于如何在地图(底图)上绘制插值站数据?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！