ngrid = 500
        x = np.array(lats)
        y = np.array(lons)
        z = np.array(data)

        xi = np.linspace(lat_min, lat_max, ngrid)
        yi = np.linspace(lon_min, lon_max, ngrid)
        xi, yi = np.meshgrid(xi, yi)
        zi = griddata(x, y, z, xi, yi)

        # draw the map
        plt.xlim(lat_min, lat_max)
        plt.ylim(lon_min, lon_max)
        plt.contour(xi, yi, zi, 20, linewidths=1)
        plt.scatter(x, y, c=z, s=20)