For concreteness sake, let's suppose we are computing the moment of inertia of the Koch Snowflake $K$ of uniform density and mass 1 shown on the left below. This snowflake is centered at the origin and has diameter 3 - that is, the maximum distance between two points is three. From here, it would be easy to compute the moment of another Koch flake of different size and uniform density, as long as it's centered at the origin.
The moment of inertia of this flake $K$, as I understand it, is simply
$$\iint\limits_{K} (x^2+y^2) \, dA$$
divided by the area of the snowflake.
There are any number of tools that can compute this numerically. Using Mathematica, I came up with 0.813719. (I could provide code, if desired.)
This value can also be computed exactly using the theory of self-similar integration as described in Bob Strichartz paper Evaluating Integrals Using Self-similarity. Applying this, I computed an exact value of $9/11$ or $0.\overline{81}$. The remainder of this answer will describe that process, though it does require a fairly sophisticated knowledge of self-similarity and measure theory.
We will use the fact that the Koch Snowflake is self-similar - composed of 7 copies of itself as shown in the figure above on the right. The scaling factor for the middle piece is $1/\sqrt{3}$ and it's $1/3$ for the others. Let us suppose that the functions in the IFS mapping $K$ onto the individual pieces are $T_0,T_1,\ldots,T_6$, with $T_0$ mapping onto the central piece. The exact functions in the IFS are
\begin{align}
T_0(x,y) &= \frac{1}{\sqrt{3}}R\left(\frac{\pi}{2}\right)\left(\begin{array}{c}{x\\y}\end{array}\right) \\
T_{i}(x,y) &= \frac{1}{3}\left(\begin{array}{c}{x\\y}\end{array}\right) +
\left(\begin{array}{c}{\cos\left(\frac{\pi}{6} + i\frac{\pi}{3}\right)\\
\sin\left(\frac{\pi}{6} + i\frac{\pi}{3}\right)}\end{array}\right),
\end{align}
where $\varphi$ ranges from $\pi/6$ to $11\pi/6$ in steps of $\pi/3$ and $R(\pi/2)$ represents a rotation through the angle $\pi/2$.
As explained by Strichartz, the integral of a function $f$ defined on a self-similar set can be computed with respect to any self-similar measure $\mu$. Specializing our notation somewhat to deal with this specific problem, a measure $\mu$ on $\mathbb R^2$ is called self-similar with respect to a list of weights $p_0,p_1,\ldots,p_6$, if it satisfies
$$\mu(A) = \sum_{i=0}^6 p_i \mu(T_i^{-1}(A)$$
for every $A\subset\mathbb R^2$. As a result, any integral with respect to $\mu$ will satisfy
$$\int f d\mu = \sum_{i=0}^6 p_i \int f \circ T_i \, d\mu.$$
As explained in Strichartz' paper there is a unique measure $\mu$ that satisfies the self-similarity condition for any given list of probabilities $p_0,p_1,\ldots,p_6$. We need to choose the probabilities so that the measure $\mu$ is uniform. Given an IFS with scaling factors $r_0, r_1, \ldots, r_6$, a uniform measure can be constructed by choosing $p_i = r_i^s$, where $s$ satisfies Moran's equation
$$\sum_{i=0}^6 r_i^s = 1.$$
This yields the probability list $p_0=1/3$ and $p_i = 1/9$ for $i=1,\ldots,6$.
Now, since the total mass of $K$ with respect to $\mu$ is 1, we certainly have that the integral of any constant is just that constant, i.e.
$$\int c \, d\mu = c.$$
We can use the self-similarity of the integral to compute the integrals of $f_1(x,y)=x$ and $f_2(x,y)=y$ as follows. Written down as a list of functions that return ordered pairs, the IFS is
\begin{align}
T_0(x,y) &= \left(-\frac{y}{\sqrt{3}},\frac{x}{\sqrt{3}}\right) \\
T_1(x,y) &= \left(\frac{x}{3}+\frac{\sqrt{3}}{2},\frac{y}{3}+\frac{1}{2}\right) \\
T_2(x,y) &= \left(\frac{x}{3},\frac{y}{3}+1\right) \\
T_3(x,y) &= \left(\frac{x}{3}-\frac{\sqrt{3}}{2},\frac{y}{3}+\frac{1}{2}\right) \\
T_4(x,y) &= \left(\frac{x}{3}-\frac{\sqrt{3}}{2},\frac{y}{3}-\frac{1}{2}\right) \\
T_5(x,y) &= \left(\frac{x}{3},\frac{y}{3}-1\right) \\
T_6(x,y) &= \left(\frac{x}{3}+\frac{\sqrt{3}}{2},\frac{y}{3}-\frac{1}{2}\right)
\end{align}
Note that $f_1$ returns just the first component and $f_2$ returns the second. Extracting those first components, multiplying them by the terms in the probability list and adding them up, to apply the self-similar integration identity, we get
$$\int x \, d\mu = \frac{2}{9}\int x \, d\mu - \frac{1}{3\sqrt{3}} \int y \, d\mu.$$
Similarly,
$$\int y \, d\mu = \frac{1}{3\sqrt{3}} \int x \, d\mu + \frac{2}{9}\int y \, d\mu.$$
This leads to a pair of equations in the unknowns
$$\int x \, d\mu \: \text{ and } \int y \, d\mu,$$
with a unique solution that both are zero. From the symmetry of those functions and the domain, this is absolutely correct. A similar procedure may be carried out on the second order terms. We find that
\begin{align}
\int x^2 \, d\mu &= \int (1/3 + 2x^2/27 +y^2/9) \, d\mu \\
\int xy \, d\mu &= -\int xy/27 \, d\mu \\
\int y^2 \, d\mu &= \int (1/3 + x^2/9 + 2y^2/27) \, d\mu.
\end{align}
Solving these for the integrals of $x^2$ and $y^2$ and adding, we obtain the desired result.
Best Answer
Let's focus first on the Koch curve. A standard construction for the curve starts with a single line segment, breaks it into thirds, and replaces the middle third with the other two sides of an equilateral triangle.
Now, by "vertex of the Koch curve", I suppose you mean one of the vertices in one of these polygonal approximations. The easiest way to prove that there are countably many of these is to simply point out that there are finitely many at each step and that the total collection of vertices is the countable union of these finite sets. Nonetheless, it is an interesting problem to enumerate them explicitly and potentially useful as there are certainly self-similar sets where these points play a distinguished role in the limit set.
The construction of the Koch curve can be described using an iterated function system $\{T_0,T_1,T_2,T_3\}$, where \begin{align} T_0(x) &= x/3 \\ T_1(x) &= R(\pi/3)\,x/3 + \langle 1/3,0 \rangle \\ T_2(x) &= R(-\pi/3)\,x/3 + \langle 1/2, \sqrt{3}/6 \rangle \\ T_3(x) &= x/3 + \langle 2/3, 0 \rangle. \\ \end{align} In this notation, $x\in\mathbb R^2$ and $$ R(\theta) = \left( \begin{array}{cc} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{array}\right) $$ is the rotation matrix through the angle $\theta$ about the origin.
Now, the segments in the first step of the approximation are exactly the images of the unit interval $((x,0):x\in [0,1])$ under the functions in the IFS. We might represent this like so:
The side labeled $\{i\}$ is the image of the unit interval under the function $T_i$. The process extends naturally to the next level.
The segment labeled $\{i_2,i_1\}$ is the image of the unit interval under the function $T_{i_1}\circ T_{i_2}$. These segments come in a natural order to form a continuous path moving from $(0,0)$ to $(1,0)$. As such, each has an initial endpoint and a terminal endpoint. The labels on the edges can be used to enumerate the vertices. Specifically, each label is the base four expansion of an integer; we'll map that integer to the initial endpoint of the corresponding segment. At level three, this produces an enumeration that looks something like so:
To map onto the vertex of the snowflake itself, you can simply map the multiples of three to the vertices above, the numbers of the form $3n+1$ to one of the other sides and the numbers of the form $3n+2$ to the remaining side. Note that the right most point above is not taken. That's fine, as that can be the initial point of the next side.