插值网格数据把数据往北转移，这是一个错误吗？

Question

我观察到来自scipy.interpolate.griddata的意外结果。我试图用matplotlib.basemap和scipy.interpolate.griddata来可视化一组不规则间隔的点。在

数据分为三个列表：纬度、经度和数值。为了让它们出现在地图上，我将数据插值到一个规则的网格上，并使用Basemap的imshow函数将其可视化。在

我观察到插值数据从真实位置向北移动。在

这里有一个例子。在这里，我要强调一个由两条经线和两条平行线组成的细胞。我希望得到这样的结果：

但是我得到的是这样的：

你可以看到红色矩形明显向北移动。在

我试图改变网格分辨率和点的数量，但是这似乎对观察到的偏移没有任何影响。在

这里有一个IPython notebook来说明这个问题。在

下面是完整的代码：

import numpy as np
from numpy import random
from scipy import interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# defining the region of interest
r = {'lon':[83.0, 95.5], 'lat':[48.5,55.5]}
# initializing Basemap
m = Basemap(projection='merc', 
            llcrnrlon=r['lon'][0],
            llcrnrlat=r['lat'][0],
            urcrnrlon=r['lon'][1],
            urcrnrlat=r['lat'][1],
            lon_0=r['lon'][0], 
            ellps='WGS84',
            fix_aspect=True,
            resolution='h')
# defining the highlighted block
block = {'lon':[89,91],'lat':[50.5,52.5]}
# generating the data
npixels = 100000
lat_range = r['lat'][1] - r['lat'][0]
lats = lat_range * random.random(npixels) + r['lat'][0]
lon_range = r['lon'][1] - r['lon'][0]
lons = lon_range * random.random(npixels) + r['lon'][0]
values = np.zeros(npixels)
for p in range(npixels):
    if block['lat'][0] < lats[p] < block['lat'][1] \
    and block['lon'][0] < lons[p] < block['lon'][1]:
        values[p] = 1.0 
# plotting the original data without interpolation
plt.figure(figsize=(5, 5))
m.drawparallels(np.arange(r['lat'][0], r['lat'][1] + 0.25, 2.0),
                    labels=[True,False,True,False])
m.drawmeridians(np.arange(r['lon'][0], r['lon'][1] + 0.25, 2.0), 
                    labels=[True,True,False,True])
m.scatter(lons,lats,c=values,latlon=True,edgecolors='none')
# interpolating on the regular grid
nx = ny = 500
mapx = np.linspace(r['lon'][0],r['lon'][1],nx)
mapy = np.linspace(r['lat'][0],r['lat'][1],ny)
mapgridx,mapgridy = np.meshgrid(mapx,mapy)
mapdata = interpolate.griddata(list(zip(lons,lats)),values,
                   (mapgridx,mapgridy),method='nearest')
# plotting the interpolated data
plt.figure(figsize=(5, 5))
m.drawparallels(np.arange(r['lat'][0], r['lat'][1] + 0.25, 2.0),
                    labels=[True,False,True,False])
m.drawmeridians(np.arange(r['lon'][0], r['lon'][1] + 0.25, 2.0), 
                    labels=[True,True,False,True])
m.imshow(mapdata)

0.0.0我在看这个

Answer 1

保利·维塔宁在SciPy bugtracker上answered the question。在

如果将basemap.imshow()替换为matplotlib.pyplot.pcolormesh()，则问题就消失了

替换上面的

m.imshow(mapdata)

与

^{pr2}$

生成正确对齐的图像。在

目前还不清楚我在basemap.imshow上做错了什么，但这可能是另一个问题。在