This is much simpler than you think. You need to understand that what underlies probability is integration. So, if you are interested in the marginal of the $k$th coordinate, what you really are looking for is knowing the value of $\int_A \pi(x_k) dx_k$ for all (measurable, whatever) $A$. How on earth would you estimate it? The solution is simple: you have a bunch of samples $X^{i} \in \mathbb{R}^d$. You estimate the above integral by counting how many times you have $X^{i}_k \in A$ and then divide by the number of samples you have - $N$. That is it.
In the notation you use, $\delta_{X_k^{i}}( \cdot)$ is a "function" s.t. if you integrate it over a set $A$, the result is 1if $X_k^i$ is in $A$ and zero otherwise: $\int_A \delta_{X_k^i}(x) dx = 1$ iff $X_k^i \in A$ and $0$ otherwise.
$$
\int_A \frac{1}{N}\sum_i \delta_{X_k^i}(x)dx = \frac{1}{N} \sum \int_A \delta_{X_k^i}(x
)dx\\
=\frac{1}{N} \cdot \text{number of times $X_k^i$ is in $A$ }
$$
just like our intuition would expect!!