How was the GCM polynomial found? As far as I understand, there is no general way to enumerate irreducible polynomials in a particular finite field, which are similar in nature to prime numbers over the integers.

The GCM mode finite field has order $2^{128}$, which matches the block size of AES, and uses the irreducible polynomial (see source page 8)

$$x^{128} + x^7 + x^2 + x + 1$$

Why does the source PDF not go into any detail regarding the choice and how was it found? It does not even give proof that the polynomial is irreducible.
For the GCM mode polynomial, it's likely that they simply looked it up in a table.  Low-weight irreducible polynomials over ${\rm GF}(2)$ are useful enough that people have spent time compiling lists of them; the one I linked to above (Seroussi 1998) is fairly often cited, and indeed contains the GCM polynomial.

Of course, this just changes the question to how such lists are compiled in the first place.  As the paper linked above notes, the basic procedure is to simply generate successive low-weight polynomials and test them for irreducibility.

The specific test that Seroussi used to generate the linked table was apparently based on Victor Shoup's NTL library, which nowadays can be found here; the GF2XFactoring module in that library provides an IterIrredTest() function that implements "an iterative deterministic irreducibility test, based on DDF."  Conveniently, there are also BuildIrred() and BuildSparseIrred() functions for finding "canonical" irreducible polynomials of a given degree over ${\rm GF}(2)$.

For more details on how irreducibility testing actually works, you might want to take a look at e.g. this answer on math.SE and these slides linked from it (Brent &amp; Zimmermann 2008).  In particular, according to the slides, it would seem that for most practical purposes (i.e. for polynomials of degree less than 100 million or so), a perfectly good way to test if a polynomial $p(x)$ of degree $d$ is irreducible over ${\rm GF}(2)$ is to simply use repeated squaring and modular reduction to test whether $$x^{2^d} \equiv x \pmod{p(x)}.$$

(Ps. There's apparently also a Wikipedia article on factorization of polynomials over finite fields that you might find useful.)