Extracting the Inner Workings of a Decision Tree I'm training an sklearn decision tree in python, and it has achieved around 99% accuracy on cross validated test data, so I would like to know what exactly the tree is doing in order for it to predict the test data so accurately.
I know I can plot feature importance and look at SHAP values, but I'm wondering whether there is any way I can look at individual nodes on the tree to find out exactly what it's doing?

You could visualize the decision tree using Graphviz

which sklearn shows how to do like this:
from sklearn.datasets import load_iris
from sklearn import tree
iris = load_iris()
X, y = iris.data, iris.target
clf = tree.DecisionTreeClassifier()
clf = clf.fit(X, y)
dot_data = tree.export_graphviz(clf, out_file=None, 
                     feature_names=iris.feature_names,  
                     class_names=iris.target_names,  
                     filled=True, rounded=True,  
                     special_characters=True)  
graph = graphviz.Source(dot_data)  
graph.view() 

Each node contains the following information:

*

*decision rule (e.g. $\text{petal length (cm)} \leq 2.45$)

*Gini impurity, $\sum_i p_i (1 - p_i)$

*number of samples that reach the node (e.g. samples=54)

*counts of members of each class (e.g. value = [0,2,4])

*dominant (or at least a non-dominated) class (e.g. class=versicolor)

And you can also visualize the decision boundary

which sklearn shows us you can do like this:
import numpy as np
import matplotlib.pyplot as plt

from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
from sklearn.inspection import DecisionBoundaryDisplay

iris = load_iris()


# Parameters
n_classes = 3
plot_colors = "ryb"
plot_step = 0.02


for pairidx, pair in enumerate([[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]]):
    # We only take the two corresponding features
    X = iris.data[:, pair]
    y = iris.target

    # Train
    clf = DecisionTreeClassifier().fit(X, y)

    # Plot the decision boundary
    ax = plt.subplot(2, 3, pairidx + 1)
    plt.tight_layout(h_pad=0.5, w_pad=0.5, pad=2.5)
    DecisionBoundaryDisplay.from_estimator(
        clf,
        X,
        cmap=plt.cm.RdYlBu,
        response_method="predict",
        ax=ax,
        xlabel=iris.feature_names[pair[0]],
        ylabel=iris.feature_names[pair[1]],
    )

    # Plot the training points
    for i, color in zip(range(n_classes), plot_colors):
        idx = np.where(y == i)
        plt.scatter(
            X[idx, 0],
            X[idx, 1],
            c=color,
            label=iris.target_names[i],
            cmap=plt.cm.RdYlBu,
            edgecolor="black",
            s=15,
        )

plt.suptitle("Decision surface of decision trees trained on pairs of features")
plt.legend(loc="lower right", borderpad=0, handletextpad=0)
_ = plt.axis("tight")

The sklearn docs uses standard orthogonal projections above. In principle you can look at other projections which may first involve rotation or other operations.
The former visualization can gives a procedural sense of what a decision tree is doing, and the latter can give you a spatial sense of how the model has partitioned the space.
