2

I am trying to interpolate sparse data over a meshgrid, but am observing some rather odd behavior. The white dots are precisely where I have values, and I am relying on the linear interpolation algorithm to fill in the other grids where possible. I recognize that this type of interpolation is not perfect due to the obvious lack of data, but how come some of the points where I have data fall outside the meshgrid that I am interpolating over? Is this a common phenomenon? This doesn't change even if I make the grid coarser.

I would appreciate some insight into why this happens, (perhaps how the linear interpolation works), or if there are any ways to fix this. See the red circles in the picture below for example: enter image description here
Data points provided for interpolation falling outside the meshgrid that is interpolated over

The following is some code on the interpolation that generated the gridded data. 

#mesh grid
xg = np.linspace(-130, -60, num=70)
yg = np.linspace(20,50,num=30)
Xg,Yg = np.meshgrid(xg,yg)

zg1 = griddata(points1, df2['tempratio'], (Xg, Yg), method = 'linear')

from mpl_toolkits.basemap import Basemap

lon_0 = xg.mean()
lat_0 = yg.mean()

m = Basemap(width=5000000, height=3500000,
           resolution='l', projection='stere',\
           lat_ts=40, lat_0=lat_0, lon_0=lon_0)

xm, ym = m(Xg, Yg)

cs = m.pcolormesh(xm,ym,zg1,shading='flat',cmap=plt.cm.Reds)
Zephyr
  • 11,891
  • 53
  • 45
  • 80
Houndbobsaw
  • 33
  • 4
  • thank you for the suggestion, I am not sure how far I can provide the data to make it a reproducible example, but I can certainly provide the code to the interpolation procedure and how this map was generated. I'm sorry I couldn't be more specific at this time. – Houndbobsaw Jun 16 '20 at 18:55

1 Answers1

3

griddata assigns values to the vertices of a grid, so 70x30 points. pcolormesh doesn't color vertices, but the rectangles in-between. There are only 69x29 rectangles formed by the given vertices. So, one row and one column of zg1 will be dropped. To counter that, an extra row and extra column can be added to the coordinates and shifting everything half a rectangle in each direction.

It still doesn't force griddata to include all given points, but goes a step towards the desired outcome. A denser grid can also help. (Choosing 'nearest' instead of 'linear' interpolation would fill the complete grid.)

Here is some code to illustrate what's happening:

import numpy as np
from scipy.interpolate import griddata
from matplotlib import pyplot as plt

def extend_range(x):
    dx = (x[1] - x[0]) / 2
    return np.append( x - dx, x[-1] + dx)

N = 10
points1 = np.vstack([np.random.randint(-130, -60, N), np.random.randint(20, 50, N)]).T
tempratio = np.random.randint(0, 20, N)

xg = np.linspace(-130, -60, num=15)
yg = np.linspace(20, 50, num=10)
Xg, Yg = np.meshgrid(xg, yg)

zg1 = griddata(points1, tempratio, (Xg, Yg), method='linear')

fig, axs = plt.subplots(ncols=2, figsize=(12, 4))
for ax in axs:
    ax.scatter(Xg, Yg, c=zg1, cmap='coolwarm', ec='g', s=80, zorder=2, label='griddata')
    ax.scatter(points1[:,0], points1[:,1], c=tempratio, cmap='coolwarm', ec='black', s=150, zorder=3, label='given data')

    if ax == axs[0]:
        ax.pcolormesh(xg, yg, zg1, shading='flat', cmap='coolwarm')
        ax.set_title('given x and y ranges')
    else:
        #todo: convert xg and yg to map coordinates
        ax.pcolormesh(extend_range(xg), extend_range(yg), zg1, shading='flat', cmap='coolwarm')
        ax.set_title('extended x and y ranges')
    ax.legend()
plt.show()

example plot

JohanC
  • 71,591
  • 8
  • 33
  • 66