I've always found the approach to describing finite element methods that focuses on the discrete linear system and works backward unnecessarily confusing. It is much clearer to go the other way, even if that involves a bit of mathematical notation in the beginning (which I'll try to keep to a minimum).
Assume that you are trying to solve an equation $A u = f$ for given $f$ and unknown $u$, where $A$ is a linear operator that maps functions (e.g., describing the displacement at every point $(x,y)$ in a domain) in a space $V$ to functions in another space (e.g., describing the applied forces). Since the function space $V$ is usually infinite-dimensional, this system cannot be solved numerically. The standard approach is therefore to replace $V$ by a finite-dimensional subspace $V_h$ and look for $u_h \in V_h$ satisfying $Au_h = f$. This is still infinite-dimensional due to the range space (which we'll assume for simplicity to be $V$ as well), so we just ask for the residual $Au_h-f\in V$ to be orthogonal to $V_h$ -- or equivalently, that $v_h^T(Au_h-f) = 0$ for every basis vector $v_h$ in $V_h$. If we now write the $u_h$ as a linear combination of these basis vectors, we are left with a linear system for the unknown coefficients in this combination. (The terms $v_i^TAu_j$ are exactly the entries of the stiffness matrix $K_{ij}$, and $v_j^Tf$ are the entries of the load vector. If $A$ is a differential operator, one usually performs integration by parts at some point, but this is not important here.)
None of this so far is specific to finite element methods, but applies to any so-called Galerkin method or method of weighted residuals.
The finite element method is characterized by a special choice of $V_h$: The computational domain is decomposed into a number of elements of the same basic shape (e.g., triangles; the process is often called triangulation), and the space $V_h$ is chosen such that restricted to each element, functions in $V_h$ are polynomials (e.g, linear in $x$ and $y$). Furthermore, the basis functions are chosen such that they are non-zero only in (the neighborhood of) one of the elements. The point of this choice is that you can build such a basis of $V_h$ fairly easily by finding a basis $\{\psi_j\}$ of the polynomial space on a single reference element (such as the triangle with vertices $(0,0)$, $(0,1)$ and $(1,0)$) and then using an affine transformation to map these basis functions to basis functions on each element in the triangulation. These $\psi_j$ are the shape functions. Usually, one requires that the local basis functions take the value $1$ at only one of the vertices and $0$ at the others (called a nodal basis), which is what the page you linked is talking about. Like any polynomial, these are uniquely determined by a number of interpolation conditions (e.g., a polynomial of degree $1$ on an interval is determined by its value at two distinct points; if the basis is a nodal basis, these values are taken at the vertices), which are also referred to as the degrees of freedom of the element. (The total degrees of freedom can be less than the sum of the degrees of all elements if some have to be fixed to ensure global properties of the approximation such as continuity.)
(Other choices of $V_h$ lead to other methods; in fact, there are spectral methods where the basis functions are chosen such that the stiffness matrix is the identity. Of course, there's no free lunch, so other parts of the procedure become more difficult with this basis. Similarly, there are other interpolation conditions than point evaluations, such as derivative evaluations or averages over elements or sides.)
