In my dynamics textbook, the introduction to the moment of inertia reads:
"Moment of inertia is the integral of the "second moment" about an axis of all the elements of mass dm which compose the body. For instance, the body's moment of inertia about the z axis is:" $$I=\int_{m} r^2 \mathrm{d}m$$
Where does this formula come from and what is the derivation for it?
Short answer
The moment of inertia is, by definition, the integral of the "second moment" about an axis of all the elements of mass $dm$ which compose the body. More generally, the definition of the inertia tensor $\mathbf J$, taking the center of mass to be the pivot point, is
$$\mathbf J \equiv \int_{\mathcal R} \left(\left(\mathbf \pi \cdot \mathbf \pi\right)\mathbf I - \mathbf \pi \otimes \mathbf \pi \right) \rho dv$$
where $\mathbf \pi = \mathbf x - \mathbf{\bar x}$ is the position vector of the element of mass $\rho dv$ relative to the center of mass, $\rho$ is the mass density of the body at $\mathbf x$, $dv$ is the volume the element of mass occupies, $\mathbf I$ is the identity tensor, and $\mathbf \pi \otimes \mathbf \pi$ is called the outer product of $\mathbf \pi$ and $\mathbf \pi$.
We define the moment of inertia and the inertia tensor these ways because their definitions appear in the expressions for angular momentum and angular kinetic energy as a natural consequence of generalizing Newtonian mechanics -- which only applies to point particles -- to the dynamic behavior of rigid bodies. In the most general case, the use of rotation tensors is required, but even in the simpler cases where only the definition you have provided for the moment of inertia is sufficient for the analysis of rigid-body motion, the idea of generalizing the mechanics of point particles to rigid bodies is the justification for this definition.
(For some basic information on tensors, feel free to check out this answer regarding the use of tensors in electromagnetism.)
Long answer
Most derivations for the moment of inertia of a body of mass $m$ use two assumptions
The following derivation is adapted and abridged from Intermediate Dynamics for Engineers: A Unified Treatment of Newton-Euler and Lagrangian Mechanics by Oliver O'Reilly. A more complete treatment of this derivation as well as rigid body dynamics in general, including an extensive discussion on the different types of tensors used to analyze rigid-body motion, can be found in Part Three of the text.
Consider a rigid body of mass $m$ occupying a region $\mathcal R$ whose center of mass, located at $\mathbf{\bar x}$, has a velocity $\mathbf{\bar v}$. Additionally, consider an arbitrary mass element $\rho dv$ located at $\mathbf x$ on the body and moving with velocity $\mathbf v$. Assumption 2.) ensures that the mass element will not experience translational motion relative to the center of mass -- it will only experience rotational motion relative to the center of mass. This allows us to express the relative velocity of the mass element as
$$\mathbf v - \mathbf{\bar v} = \mathbf \omega \times \left(\mathbf x - \mathbf{\bar x}\right) = \mathbf \omega \times \mathbf \pi \tag{1}$$
where $\mathbf \omega$ is the angular velocity vector of the entire body -- not just the mass element $\rho dv$ -- relative to the center of mass.
From Newtonian mechanics, the angular momentum of a mass element about the center of mass, $d \mathbf H$, by definition is
$$d \mathbf H = \mathbf \pi \times \rho dv \mathbf v = \mathbf \pi \times \rho dv \left(\mathbf{\bar v} + \mathbf \omega \times \mathbf \pi \right)\tag{2}$$
Integrating $d \mathbf H$ over the entire region $\mathcal R$ gives us the total angular momentum of the body about the center of mass, $\mathbf H$,
$$\begin{align}\mathbf H & = \int_{\mathcal R} \left(\mathbf \pi \times \left(\mathbf{\bar v} + \mathbf \omega \times \mathbf \pi \right)\right) \rho dv \\ & = \int_{\mathcal R} \left(\mathbf \pi \times \mathbf{\bar v} + \mathbf \pi \times \left(\mathbf \omega \times \mathbf \pi \right)\right) \rho dv \\ & = \int_{\mathcal R} \mathbf \pi \times \mathbf{\bar v} \rho dv + \int_{\mathcal R} \mathbf \pi \times \left(\mathbf \omega \times \mathbf \pi \right) \rho dv\end{align} \tag{3}$$
Let's consider the first term on the right-hand side of (3). Since $\mathbf{\bar x}$ and $\mathbf{\bar v}$ are independent of the region integrated, we can take them outside of the integral,
$$\begin{align} \int_{\mathcal R} \mathbf \pi \times \mathbf{\bar v} \rho dv & = \int_{\mathcal R} \mathbf \pi \rho dv \times \mathbf{\bar v} \\ & = \int_{\mathcal R} \left(\mathbf x - \mathbf{\bar x}\right) \rho dv \times \mathbf{\bar v} \\ & = \left(\int_{\mathcal R} \mathbf x \rho dv - \int_{\mathcal R} \mathbf{\bar x} \rho dv\right) \times \mathbf{\bar v} \\ & = \left(\int_{\mathcal R} \mathbf x \rho dv - \mathbf{\bar x} \int_{\mathcal R} \rho dv\right) \times \mathbf{\bar v} \\ & = \left(\int_{\mathcal R} \mathbf x \rho dv - m \mathbf{\bar x}\right) \times \mathbf{\bar v} \end{align} \tag{4}$$
But the first term on the right-hand side of (4), by the definition of center of mass, is $m \mathbf{\bar x}$, and so,
$$\int_{\mathcal R} \mathbf \pi \times \mathbf{\bar v} \rho dv = \left(m \mathbf{\bar x} - m \mathbf{\bar x}\right) \times \mathbf{\bar v} = \mathbf 0 \times \mathbf{\bar v} = \mathbf 0 \tag{5}$$
Thus, from (3),
$$\begin{align}\mathbf H & = \int_{\mathcal R} \mathbf \pi \times \left(\mathbf \omega \times \mathbf \pi \right) \rho dv \\ & = \int_{\mathcal R} \left(\left(\mathbf \pi \cdot \mathbf \pi \right) \mathbf \omega - \left(\mathbf \pi \cdot \mathbf \omega \right) \mathbf \pi \right) \rho dv \end{align} \tag{6}$$
where we have used the vector product identity $\mathbf a \times (\mathbf b \times \mathbf c) = (\mathbf a \cdot \mathbf c) \mathbf b - (\mathbf a \cdot \mathbf b) \mathbf c$.
Defining the outer product between two vectors $\mathbf a$ and $\mathbf b$, denoted $\mathbf a \otimes \mathbf b$, as the operation with the property,
$$(\mathbf a \otimes \mathbf b) \cdot \mathbf c = (\mathbf c \cdot \mathbf b) \mathbf a$$
$$\mathbf c \cdot (\mathbf a \otimes \mathbf b) = (\mathbf c \cdot \mathbf a) \mathbf b \tag{7}$$
and the identity tensor $\mathbf I$ as the second-order tensor with the property that for any vector $\mathbf a$,
$$\mathbf I \mathbf a = \mathbf a \tag{8}$$
we can rewrite (6) as,
$$\begin{align}\mathbf H & = \int_{\mathcal R} \left(\left(\mathbf \pi \cdot \mathbf \pi \right) \mathbf I \mathbf \omega - \left(\mathbf \pi \otimes \mathbf \pi \right) \mathbf \omega \right) \rho dv \\ & = \int_{\mathcal R} \left(\left(\mathbf \pi \cdot \mathbf \pi \right) \mathbf I - \mathbf \pi \otimes \mathbf \pi \right) \mathbf \omega \rho dv \\ & = \left(\int_{\mathcal R} \left(\left(\mathbf \pi \cdot \mathbf \pi \right) \mathbf I - \mathbf \pi \otimes \mathbf \pi \right) \rho dv \right) \mathbf \omega \\ & = \mathbf J \mathbf \omega\end{align} \tag{9}$$
where we have defined $\mathbf J$ as given at the top of this answer. Similar derivations will yield the expected results for the angular kinetic energy of a rigid body, the value of the inertia tensor $\mathbf J_{P}$ about any arbitrary pivot point $P$, and the parallel axis theorem.
To complement the very nice answer by JM1, the "zeroth moment"
$$ M = \int_m dm $$
is just the total mass of the system, and the "first moment"
$$ M \mathbf R_\mathrm{c.o.m.} = \int_m \mathbf r\,dm $$
gives the location of the system's center of mass.
The whole mass moment of inertia tensor comes out naturally when considering the angular momentum of a rotating rigid body.
The goal of defining a mass moment of inertia is to separate the motion and the geometry/mass parts out of angular momentum. Then the time derivative of angular momentum can be calculated in order to establish the rotational laws of motion.
In short, in a collection of $n$ particles, all rotating together, and each with mass ${\rm d}m$, and location $\mathbf{r}$. The total angular momentum is
$$ \mathbf{L} = \int \mathbf{r} \times {\rm d} \mathbf{p}$$ where ${\rm d} \mathbf{p} = \mathbf{v}\, {\rm d}m$ is the linear momentum of each particle. Since the particle is rotating then, $\mathbf{v} = {\boldsymbol \omega} \times \mathbf{r} = - \mathbf{r} \times {\boldsymbol \omega}$. So all together we have
$$ \mathbf{L} = \int (-\mathbf{r} \times (\mathbf{r} \times {\boldsymbol \omega})) \;{\rm d}m $$
Now mass moment of inertia $\mathrm{I}$ is extracted by
$$ \mathbf{L} = \mathrm{I}\, {\boldsymbol \omega} $$ $$ \mathrm{I} = \int \left| \matrix{ y^2+z^2 & -x y & - x z \\ - x y & x^2+z^2 & - y z \\ -x z & - y z & x^2+y^2} \right| \,{\rm d}m$$
In 2D, with $z=0$ and $r^2 =x^2+y^2$ the above simplifies to $$ \mathrm{I}_{zz} = \int r^2 {\rm d}m$$
The shorthand notation for the mass moment of inertia tensor is
$$ \mathrm{I} = \int -[\mathbf{r}\times][\mathbf{r}\times]\,{\rm d}m $$
where each $[\mathbf{r}\times]$ is a 3×3 skew-symmetric matrix representing the cross product operator
$$[\mathbf{r}\times] \equiv \left| \matrix{0 & -z & y \\ z & 0 & -x \\ -y & x & 0} \right|$$
