Type error when kriging with python: TypeError - Input z must be at least a (2, 2) shaped array, but has shape (1, 1)

Question

I am new to kriging with python, and I have a csv file that I am trying to interpolate and plot. However, I am getting the following error:

TypeError - Input z must be at least a (2, 2) shaped array, but has shape (1, 1)

What can I do to correct this? I'm not sure what is going on.

In addition, if anyone has any tips on how to make this code fit my data better, I would greatly appreciate it. This code was obtained from this tutorial.

Here is my code:

#Import modules
import pandas as pd
import numpy as np
import os
import glob
from pykrige.ok import OrdinaryKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap ####You may need to pip install basemap-data-hires if you cannot get this package to load
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.patches import Path, PathPatch


#Set workplace directory
os.chdir('workplace_where_file_is')


nut = pd.read_csv(csv.file)
nut.head(55)

Output:

    Sample_ID   longitude   latitude    var1    var2    var3
0   1   -100.922987 39.248554   86  5   72
1   2   -100.921920 39.248792   90  38  3
2   3   -100.920808 39.248584   15  57  37
3   4   -100.919779 39.248591   11  38  41
4   5   -100.918698 39.247933   29  24  89
5   6   -100.919457 39.247719   86  15  95
6   7   -100.920822 39.247916   49  25  75
7   8   -100.922241 39.247906   99  59  84
8   9   -100.923254 39.247781   11  78  100
9   10  -100.924217 39.247858   50  18  41
10  11  -100.925033 39.246966   21  73  51
11  12  -100.924190 39.246857   52  13  62
12  13  -100.923128 39.247089   30  34  80
13  14  -100.921795 39.246824   61  15  11
14  15  -100.920870 39.247048   57  52  62
15  16  -100.919414 39.246841   10  7   68
16  17  -100.918580 39.247120   86  51  43
17  18  -100.917600 39.246849   22  43  26
18  19  -100.917600 39.246111   10  49  50
19  20  -100.918410 39.245921   82  51  52
20  21  -100.919478 39.246194   72  46  4
21  22  -100.920773 39.245991   80  56  66
22  23  -100.921784 39.246063   97  16  30
23  24  -100.923103 39.246296   43  100 98
24  25  -100.924202 39.246112   26  54  15
25  26  -100.925255 39.246100   53  75  95
26  27  -100.925431 39.245284   13  64  46
27  28  -100.924429 39.245093   59  93  21
28  29  -100.923143 39.245226   24  60  64
29  30  -100.922008 39.245047   59  93  37
30  31  -100.920764 39.245258   63  33  92
31  32  -100.919717 39.244983   22  16  26
32  33  -100.918345 39.245352   32  98  6
33  34  -100.917390 39.245077   45  65  42
34  35  -100.917553 39.244407   77  89  44
35  36  -100.918440 39.244260   4   27  61
36  37  -100.919476 39.244483   46  63  83
37  38  -100.920955 39.244112   85  76  43
38  39  -100.921705 39.244254   53  87  79
39  40  -100.923370 39.244238   9   94  69
40  41  -100.924262 39.244383   47  88  94
41  42  -100.925087 39.244483   63  35  20
42  43  -100.924365 39.243440   4   23  67
43  44  -100.923045 39.243196   84  65  99
44  45  -100.921845 39.243471   93  1   24
45  46  -100.920755 39.243252   29  2   24
46  47  -100.919772 39.243325   6   4   33
47  48  -100.918568 39.243501   95  79  65
48  49  -100.919691 39.242702   8   75  4
49  50  -100.920656 39.242624   10  51  40
50  51  -100.922180 39.242513   33  46  68
51  52  -100.922826 39.242723   51  57  94

Generate interpolation maps for Var1

lon = np.array(nut['longitude'])
lat = np.array(nut['latitude'])
nit = np.array(nut['var1'])

grid_space = 0.01
grid_lon = np.arange(np.amin(lon), np.amax(lon), grid_space)
grid_lat = np.arange(np.amin(lat), np.amax(lat), grid_space)

#Ordinary kriging
krig = OrdinaryKriging(lon, lat, nit, variogram_model = 'gaussian', verbose = True, enable_plotting = False, nlags = 20)
z1, ss1 = krig.execute('grid', grid_lon, grid_lat)

#Plotting the interpolation
xintrp, yintrp = np.meshgrid(grid_lon, grid_lat)
fig, ax = plt.subplots(figsize = (10, 10))
m = Basemap(llcrnrlon = lon.min()-0.1, llcrnrlat = lat.min()-0.1, urcrnrlon = lon.max()+0.1, urcrnrlat = lat.max()+0.1, projection = 'merc', resolution = 'h', area_thresh = 1000., ax = ax)

m.drawcoastlines()
x,y = m(xintrp, yintrp)
ln,lt = m(lon, lat)
cs = ax.contourf(x, y, z1, np.linspace(0, 55, 1), extend = 'both', cmap = 'jet')
cbar = m.colorbar(cs, location = 'right', pad = "7%")

#Draw parallels
parallels = np.arrange(21.5, 26.0, 0.5)
m.drawparallels(parallels, labels = [1, 0, 0, 0], fontsize = 14, linewidth = 0.0)

#Draw meridians
meridians = np.arrange(119.5, 122.5, 0.5)
m.drawmeridians(meridians, labels = [0, 0, 0, 1], fontsize = 14, linewidth = 0)

The following error is produced:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [32], in <module>
      6 x,y = m(xintrp, yintrp)
      7 ln,lt = m(lon, lat)
----> 8 cs = ax.contourf(x, y, z1, np.linspace(0, 55, 1), extend = 'both', cmap = 'jet')
      9 cbar = m.colorbar(cs, location = 'right', pad = "7%")
     11 #Draw parallels

File ~\anaconda3\envs\lib\site-packages\matplotlib\__init__.py:1412, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1409 @functools.wraps(func)
   1410 def inner(ax, *args, data=None, **kwargs):
   1411     if data is None:
-> 1412         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1414     bound = new_sig.bind(ax, *args, **kwargs)
   1415     auto_label = (bound.arguments.get(label_namer)
   1416                   or bound.kwargs.get(label_namer))

File ~\anaconda3\envs\lib\site-packages\matplotlib\axes\_axes.py:6317, in Axes.contourf(self, *args, **kwargs)
   6308 """
   6309 Plot filled contours.
   6310 
   (...)
   6314 %(contour_doc)s
   6315 """
   6316 kwargs['filled'] = True
-> 6317 contours = mcontour.QuadContourSet(self, *args, **kwargs)
   6318 self._request_autoscale_view()
   6319 return contours

File ~\anaconda3\envs\lib\site-packages\matplotlib\contour.py:812, in ContourSet.__init__(self, ax, levels, filled, linewidths, linestyles, hatches, alpha, origin, extent, cmap, colors, norm, vmin, vmax, extend, antialiased, nchunk, locator, transform, *args, **kwargs)
    808     self.origin = mpl.rcParams['image.origin']
    810 self._transform = transform
--> 812 kwargs = self._process_args(*args, **kwargs)
    813 self._process_levels()
    815 self._extend_min = self.extend in ['min', 'both']

File ~\anaconda3\envs\lib\site-packages\matplotlib\contour.py:1441, in QuadContourSet._process_args(self, corner_mask, *args, **kwargs)
   1438     corner_mask = mpl.rcParams['contour.corner_mask']
   1439 self._corner_mask = corner_mask
-> 1441 x, y, z = self._contour_args(args, kwargs)
   1443 _mask = ma.getmask(z)
   1444 if _mask is ma.nomask or not _mask.any():

File ~\anaconda3\envs\lib\site-packages\matplotlib\contour.py:1480, in QuadContourSet._contour_args(self, args, kwargs)
   1478     args = args[1:]
   1479 elif Nargs <= 4:
-> 1480     x, y, z = self._check_xyz(args[:3], kwargs)
   1481     args = args[3:]
   1482 else:

File ~\anaconda3\envs\lib\site-packages\matplotlib\contour.py:1510, in QuadContourSet._check_xyz(self, args, kwargs)
   1508     raise TypeError(f"Input z must be 2D, not {z.ndim}D")
   1509 if z.shape[0] < 2 or z.shape[1] < 2:
-> 1510     raise TypeError(f"Input z must be at least a (2, 2) shaped array, "
   1511                     f"but has shape {z.shape}")
   1512 Ny, Nx = z.shape
   1514 if x.ndim != y.ndim:

TypeError: Input z must be at least a (2, 2) shaped array, but has shape (1, 1)